Rutting Prediction of Asphalt Mixtures Modified by Polypropylene Fibers via Repeated Creep Testing by Utilising Genetic Programming

A novel application of genetic programming (GP) for modelling and presenting closed form solutions to the rutting prediction for polypropylene (PP) modified asphalt mixtures is investigated. Various PP fibers have been utilised for bitumen modification and repeated creep (RC) tests have been carried out. Marshall specimens, fabricated with multifilament 3 mm (M-03) type PP fibers at optimum bitumen content of 5% have been tested under different load values and patterns at 50 °C to investigate their rutting potential. It has been shown that the service lives of PP fiber-reinforced Marshall specimens are respectively longer than the control specimens under the same testing conditions (5 to 12 times). Input variables in the developed GP model use the physical properties of Marshall specimens such as PP type, specimen height, unit weight, voids in mineral aggregate, voids filled with asphalt, air voids, rest period and pulse counts. The performance of the accuracy of the proposed GP model is observed to be quite satisfactory. To obtain the main effects plot, detailed parametric studies have been performed. The presented closed form solution will also help further researchers willing to perform studies on the prediction of the rutting potential of asphalt without carrying out destructive tests for similar type of aggregate sources, bitumen, aggregate gradation, modification technique and laboratory conditions.


Introduction
In this study, a genetic programming (GP) model has been proposed to predict the strain accumulation (rutting potential) developed in the polypropylene (PP) fiber modified Marshall specimens during repeated load creep tests.The proposed study differs from the previous researches in the sense that, no studies have been carried out for the prediction of the RC test results with the aid of GP, closed form solutions and parametric studies carried out on Marshall specimens fabricated in the laboratory environment up to date.
The first part of this study gives short information on creep (basically repeated) testing of asphalt spanning in the last two decades about the actual loading simulation efforts.Then, well known applications of GP techniques in pavement engineering are explored.Next, experimental program presenting the results of the RC tests has been stated in a detailed manner.At this point, a very general overview of GP is presented.Expression tree for the strain accumulation under various RC loading patterns is given in the next section.In order to obtain the main effects plot, a wide range of parametric studies have been performed by using the GP models.The analysis of these results is presented next.Finally, the MATLAB code for the closed form solutions of the prediction of rutting potential is given in Appendix 1.

The Idea Behind the Repeated Creep
Testing Procedure and Recent Literature In the undertaken RC studies, first of all, the testing temperature was chosen as 50 °C to simulate actual in-situ conditions.An axial stress, σ, of 500 kPa was applied to the specimens until the specimen enters the tertiary creep region to nearly a failure point to simulate the actual inplace conditions in a realistic manner [1][2][3] .Universal testing system (UTM) was utilised in order to log the accumulation mechanisms of the developing strains or in other words rutting potential.Prior to testing, the specimens were put into the chamber of UTM for 24 hours in order to have uniform temperature distribution.To understand the behaviour of the asphalt specimens under different loading patterns, different constant stress values were chosen (namely 100, 207 and 500 kPa).As PP fiber modification was carried out, utilizing lower stress values like 100 and 207 kPa was not feasible, since under such loading the tertiary creep region could not be observed within a reasonable period of time 3 .Therefore, in order to be able to differentiate between the reference and fiber-reinforced samples, a real destructive loading level of 500 kPa (approximately 73 psi) was chosen as the standard stress value which is a main departure from the published pioneering literature of the rule of thumbs about creep testing 3 .This value very well represents the actual tire pressure of a loaded truck.The specimen strain during the pulsed loading stage of the test were measured in the same axis as the applied stress using two linear variable displacement transducers (LVDTs).The applied force was open loop controlled and rectangular in shape.Load periods were chosen as 500 ms for all of the specimens and the rest periods were 500, 1000, 1500 and 2000 ms, respectively.Four specimens were tested for each loading pattern 3 .
Matthews and Monismith 4 have performed unconfined creep tests at temperatures 25 °C, 38 °C and 49 °C which is a main departure from the published literature up to date in the testing temperature manner.In another study by Mallick et al. 5 , in order to simulate the average pavement temperature throughout the United States, 60 °C of testing temperature was utilised.Ramsamooj and Ramadan 6 had carried out creep tests at four stress levels under constant stresses of 150, 400, 650 and 900 kPa.Zhang et al. 7 had utilised a material testing system to conduct RC tests.The test temperature was 60 °C.Test loading consisted of a 138 kPa (20 psi) confining pressure and an 827 kPa (120 psi) normal pressure.Tashman et al. 8 had carried out triaxial confined static creep test in determining the model parameters related to their studies.Chen et al. 9 had investigated the mechanical responses and modelling of rutting in flexible pavements.Goh and You 10 has revised the RC testing ambient temperature in a different manner than only utilising 40°C.Vardanega et al. 11 have used a very similar ambient temperature and loading pattern to Tapkın et al. in the studies that they have carried out (namely 50 °C and 0.5 (also 1.5) second loading and 1.5 (also 0.4) second unloading) depending on the argument which is also verified by the experience of QDMR (Queensland Department of Main Roads).Chen et al. 12 investigated the utilization of recycled brick powder as alternative filler in asphalt mixture.They had carried out static and dynamic creep tests using UTM.

Recent Genetic Programming Studies Carried Out in Pavement Engineering Applications
Sundin and Braban-Ledoux 13 have summarized the findings of up to date research articles concerning the application of artificial intelligence to pavement management and to illustrate the potential of such tools.Hadi and Arfiadi 14 have mentioned about the design of rigid pavements according to AUSTRDADS.A genetic algorithm is used to find the optimum design.Chan et al. 15 have tried to solve the problem of pavement maintenance management at the network level and genetic algorithms (GAs) have been demonstrated to be better than traditional techniques in terms of solution quality and diversity.Tack and Chou 16 have described in their studies, the development of a GA based optimization tool for determining the optimal multiyear pavement repair schedule.Tsai et al. 17 have demonstrated the applicability of the GA to solve nonlinear optimization problems encountered in asphalt pavement design.Tsai et al. 18 , furthermore, inspected the three stage Weibull equation and tree based model to characterize the mix fatigue damage process.In another study by Alavi et al. 19 , a high precision model was derived to predict the flow number of dense asphalt mixtures using a novel hybrid method coupling GP and simulated annealing.Gandomi et al. 20 also developed a novel hybrid method coupling GP and orthogonal least squares.In another paper by Gandomi and Alavi 21 , a new approach for behavioural modelling of structural engineering systems is presented using a promising variant of GP, namely multigene GP (MGGP).In the second part of their studies, this time, Gandomi and Alavi 22 present an endeavour to exploit a robust MGGP method for the analysis of geotechnical and earthquake engineering systems.Finally, Jorge and Ferreira 23 present a new maintenance optimisation system, called GENEPAV-HDM4, which was developed to integrate the pavement management system of the municipality of Viseu (Portugal).

Experiments Carried Out
Continuous aggregate gradation has been used to fit the gradation limits for wearing course Type 2 set by General Directorate of Turkish Highways 24 .The aggregate was a calcareous type crushed stone obtained from a local quarry and 50/70 penetration bitumen was obtained from a local refinery were used.Physical properties of the bitumen samples are given in Table 1.The physical properties of coarse and fine aggregates are given in Tables 2 and 3.The apparent specific gravity of filler is 2790 kg.m -3 .The mixture gradation and gradation limits are given in Table 4.In the wet basis modification procedure, standard 50/70 penetration bitumen was modified by utilising PP fibers.The fibers were premixed with bitumen using a standard mixer at 500 revolutions per minute (rpm) for two hours.This low shear rate is one of the main advantages of PP fiber modification when compared to other types of polymeric bitumen modifiers which need high shear mixing rates to attain the required the compatibility between the bituminous binder and the modifier.The mixing temperature was around 165-170 °C25 .According to the workability criteria, multifilament 3 mm (M-03) type fibers were found to be the best modifiers and due to the consistency of the Marshall test results, 3‰ fiber content was determined as the optimal addition amount by weight of aggregate.With these amounts, PP fibers melt in bitumen and bitumen forms a continuous phase.The physical properties of the PP fiber based bitumen samples with 3‰ fiber content by weight are given in Table 5.
The performance characteristics, such as penetration, penetration index, ductility, loss on heating, specific gravity, and softening point of the fiber modified bitumen samples were greatly improved as compared to reference specimens given in Table 1.Therefore, the addition of 3‰ of M-03 type fibers clearly shows the decrease in temperature susceptibility of the reference bitumen (as shown by the eminent increase in the penetration index of PP modified bitumen samples) providing the most significant effect on the properties of resultant asphalt concrete mixtures as an increase in the stiffness values.
To determine the optimum bitumen content, the bitumen contents corresponding to the mixtures with maximal stability and unit weight, 4% air voids and 70% voids filled with asphalt, were found and averaged 24 .These optimum bitumen contents are represented in Figure 1.
Based on these results, M-03 PP fibers at a dosage of 3‰ by the weight of aggregates were selected as optimal PP addition amount.Also, it can be seen that the optimum bitumen contents for reference and specimens with 3‰ of M-03 fibers are 4.81% and 4.97%, respectively (Figure 1).For the next step of experiments, these two values were taken as 5% for the sake of ease in preparation of the reference and modified asphalt specimens 1 .

Repeated creep tests performed
The reference specimens utilised in the RC testing were prepared with 5% optimum bitumen content.The fiber-reinforced (M-03 type with dosage of 3‰ by the weight of aggregate) specimens were also prepared with 5% bitumen content.
The RC test results are visualized through Figures 2-5.These sets of graphs present the log of strain accumulations developed versus pulse counts.These graphs show the general trend of the four similar (in lieu of proximity of air void values concept) specimens under the specified loading and temperature conditions 3 .For the sake of uniformity reasons, the constant stress had been taken as 500 kPa and the test temperature was 50 °C for all of the experiments As can be seen from Figures 2 to 5, the service lives of fiber-reinforced specimens are respectively longer than the control specimens (5 to 12 times).This is a very significant difference showing the positive effect of PP fiber modification.As can be seen from all figures, the control specimens are entering to the tertiary stage of creep only at around 2000 pulse counts.This loading rate corresponds to the primary creep stage for the PP fiber modified specimens.Fiber modified specimens reach the tertiary creep stage at much higher pulse counts.At the end of the RC tests, the control specimens have a total collapse, while the fiber modified specimens did not show any sign of failure (so, the fiber modified specimens would have had even a longer service life if the tests were continued).Therefore it can be easily  concluded that PP modification substantially increases the service lives of asphalt specimens under repeated loading.Rutting problems are not encountered in the PP modified specimens until very late pulse counts.This is valid for all of the loading patterns covered in this study.But the question is: can these strain accumulation values be estimated by some means so that the rutting potential of asphalt mixes can be predicted in an accurate manner?This study tries to find a solution to this problem from this point forward.

Background in Genetic Programming (GP)
Genetic algorithm (GA) is an optimization and search technique which is mainly based on the principles of genetics and natural selection.A GA allows a population composed of many individuals to evolve under specified selection rules to a state that maximizes the "fitness".The method was developed by John Holland and finally popularized by one of his students 26,27 .The fitness of each individual in a GA is the measure of the individual that has been adapted to the problem which is solved by employing this individual.The basis of genetic algorithms is the selection of individuals in accordance with their fitness; thus, fitness is obviously a critical criterion for optimization 28 .
GP is an extension of GAs proposed by Koza 29 .Koza defines GP as a domain independent problem solving approach in which computer programs are evolved to solve, or approximately solve the problems based on the Darwinian principle of reproduction and survival of the  fittest and analogs of naturally occurring genetic operations such as crossover (sexual recombination) and mutation.GP reproduces computer programs to solve problems whose flowchart can be found in the relevant literature 29 .

Brief overview of gene expression programming (GEP)
Gene expression programming (GEP) software (which is utilised in this study) is an extension of GP that evolves computer programs of different sizes and shapes encoded in linear chromosomes of fixed length.The chromosomes are composed of multiple genes.Each gene is encoding a smaller sub program.Furthermore, the structural and functional organization of the linear chromosomes allows the unconstrained operation of important genetic operators such as mutation, transposition, and recombination.APS 3.0, a GEP software developed by Candida Ferreira has been utilised in this study 30 .
Thus, the two main parameters of GEP are the chromosomes and expression trees (ETs).The process of information decoding (from the chromosomes to the ETs) is called translation which is based on a set of rules.The genetic code is very simple where there exist one to one relationships between the symbols of the chromosome and the functions or terminals they represent.The rules, which are also very simple, determine the spatial organization of the functions and terminals in the ETs and the type of interaction between sub ETs [31][32][33] .That is why two languages are utilized in GEP: the language of the genes and the language of ETs.A significant advantage of GEP is that it enables to infer exactly the phenotype given the sequence of a gene, and vice versa which is termed as Karva language.Consider, for example, the algebraic expression (d4* ( ) -d4) that can be represented by a diagram (Figure 6) which is called the expression tree.
For each problem, the type of linking function, as well as the number of genes and the length of each gene, is a priori chosen for each problem.While attempting to solve a problem, one can always start by using a single gene chromosome and then proceed with increasing the length of the head.If it becomes very large, one can increase the number of genes and obviously choose a function to link the sub ETs.Dne can start with addition for algebraic expressions or for Boolean expressions, but in some cases another linking function might be more appropriate (like multiplication or IF, for instance).The idea, of course, is to find a good solution, and GEP provides the means of finding one very efficiently 33 .Some detailed illustrative examples considering showing how GEP can be used to model complex realities can be found in the relevant literature [34][35][36] .

Solving a simple problem with GEP
As an illustrative example, consider the following case where the objective is to show how GEP can be used to model complex realities with high accuracy.So, suppose one is given a sampling of the numerical values from the curve (remember, however, that in real-world problems of the functions are obviously unknown): Dver 10 randomly chosen points in the real interval [-10, +10], the aim is to find a function fitting those values within a certain error.In this case, a sample of data in the form of 10 pairs (a i , y i ) is given where a i is the value of the independent variable in the given interval and y i is the respective value of the dependent variable (a i values: -4.2605, -2.0437, -9.8317, -8.6491, 0.7328, -3.6101, 2.7429, -1.8999, -4.8852, 7.3998; the corresponding y i values can be easily evaluated).These 10 pairs are the fitness cases (the input) that will be used as the adaptation environment.The fitness of a particular program will depend on how well it performs in this environment 31 .
There are five major steps for preparing to use gene expression programming.The first is to choose the fitness function.For this problem one could measure the fitness f i of an individual program i by the following expression: where M is the range of selection, C (i,j) the value returned by the individual chromosome i for fitness case j (out of C t fitness cases) and T j is the target value for fitness case j.If, for all j, ( ) T (the precision) is less than or equal to 0.01, then the precision is equal to zero, and For this problem, M = 100, is used therefore, f max = 1000.The advantage of this kind of fitness function is that the system can find the optimal solution for itself.However there are other fitness functions available which can be appropriate for different problem types 31 .
The second step is choosing the set of terminals T and the set of functions F to create the chromosomes.In this problem, the terminal set consists obviously of the independent variable, i.e., T = {a}.The choice of the appropriate function set is not so obvious, but a good guess can always be made in order to include all the necessary functions.In this case, to make things simple, the four basic arithmetic operators are used.Thus; F= {+, -, *, /}.It should be noted that there many other functions that can also be utilised.
The third step is to choose the chromosomal architecture, i.e., the length of the head and the number of genes.
The fourth major step for preparing to use gene expression programming is to choose the linking function.In this case we will link the sub ETs by addition.Dther linking functions are also available such as subtraction, multiplication and division.
Finally, the fifth step is to choose the set of genetic operators that cause variation and their rates.In this case one can use a combination of all genetic operators (mutation at p m = 0.051; IS and RIS transposition at rates of 0.1 and three transposons of length 1, 2, and 3; one-point and two-point recombination at rates of 0.3; gene transposition and gene recombination both at rates of 0.1).
To solve this problem, let's choose an evolutionary time of 50 generations and a small population of 20 individuals in order to simplify the analysis of the evolutionary process and not to fill up this manuscript with pages of encoded individuals.However, one of the advantages of GEP is that it is capable of solving relatively complex problems using small population sizes and, thanks to the compact Karva notation that it is possible to fully analyse the evolutionary history of a run.
A perfect solution can be found in generation 3 which has the maximum value 1000 of fitness.The sub ETs codified by each gene are given in Figure 7.Note that it corresponds exactly to the same test function given above in Equation 1 31 .
Thus expressions for each corresponding Sub-ET can be given as follows: The main focus of this study is to explore an application of GP for modelling and presenting closed form solutions to the rutting prediction with the aid of RC testing results for PP modified asphalt.Therefore an extensive laboratory testing and data analysis phase has been performed.The details of the experimental database including the ranges of parameters are given in Table 6.Prior to GP modelling, various regression models had been obtained in order to be able to observe the correlation between input and output variables.Therefore multivariate linear, nonlinear and stepwise regression models were obtained which are shown in Table 7.As seen, the performances of the models are extremely poor when compared to GP modelling.Regarding nonlinear regression modelling, only two models could be obtained for three and four variables.Dther more complex nonlinear models with more than four inputs did not lead to any significant result.

Results of GP formulations
The input variables in the developed GP model use the physical properties of standard Marshall specimens such as PP type, specimen height, unit weight, voids in mineral aggregate, voids filled with asphalt, air voids and RC test properties such as rest period and pulse counts in order to predict the rutting potential of the fabricated specimens.Prior to GP modelling, the experimental results are divided into randomly selected training and testing sets among the experimental database as 80% and 20% respectively to prevent over fitting.The GP modelling was performed by GeneXpro-Tools 4 30 .The GP model was constructed with training sets and the accuracy was verified by testing sets which the GP model is facing up with for the first time.Related parameters for the training of the GP models are given in Table 8.It should be noted that the proposed GP formulation is valid for the ranges of training set given in Table 6 and similar type of aggregate sources, bitumen, aggregate gradation, mix proportioning, modification technique and laboratory conditions.Statistical parameters of test and training sets of GP formulations are presented in Table 9.The GP results versus actual test results for the strain accumulation are represented in Figure 8.

Parametric Studies
The main effects plot is an important graphical tool to visualize the independent impact of each variable utilised on strain accumulation values.This tool allows the reader to visualise a better and much simpler snapshot of the overall significance of variable effects on the outputs.In the main effects plot, the mean output is plotted at each factor level which is later connected by a straight line.The slope of the line for each variable is the degree of its effect on the output.In order to obtain the main effects plot, a wide range of detailed parametric studies have been performed by utilising the proposed GP model.The mean values of all variables are used to observe the general trend of each variable.The evaluation of separate interaction effect plots between any two variables is also presented by using the mean values of all variables.The main effects plot will also help further researchers willing to perform studies on rutting potential without carrying out destructive tests for similar type of aggregate sources, bitumen, aggregate gradation, mix proportioning, modification technique and laboratory conditions.

Analysis of results
In this part of the study, analyses of the graphs that have been obtained at the end of parametric studies have been given.
In Figure 10 and the rest of the figures, the y-axes represent the strain accumulation values in µε.In the above figure, it can be visualized that as the pulse counts increase, the strains developing in the Marshall specimen bodies also increase in a considerable manner.This fact is also verified by the four different types of loading pattern (namely 500 ms load 500 ms rest, 500 ms load 1000 seconds rest, 500 ms load 1500 ms rest and 500 ms load 2000 ms rest).
Figure 11 presents the same phenomenon from another approach.This time, one can easily notify that at smaller number of pulse counts, strain accumulation is not affected in a considerable manner.By the end of the test period, the strain accumulation values reach to 50000 µε and more.This is an expected case from the materials engineering point of view.
As air voids start to increase in an abrupt manner, the strain accumulation obeys this trend.The reader can visualise that the change in the strain accumulation is only in a small band of 1200 µε.The change in the specimen heights does not alter the whole picture (Please refer to Figure 12).Similar arguments are valid for Figure 13.The increase in the V f (voids filled with asphalt) values mean that the specimens have lesser total air voids (Please refer to Figure 14).By moving from this point on, the decrease in the strain accumulation can be explained easily.At the first glance, it might seem that the unit weight change is not coupled with the above discussion but one must not forget that the strain accumulation band is only 70 µε.Figure 15 offers another approach to the same analysis results.
In Figure 16, as the loading pattern becomes more aggressive, strain accumulation increases in a noticeable manner.
The increase in V.M.A. (voids in mineral aggregate) values coupled with the more aggressive loading patterns result in greater strain accumulation (Please refer to Figure 17).Dne has to bear in mind that the 16% V.M.A. value is a sort of inflection point for the eminent increase of strain accumulation which stands for the tertiary creep region.Finally Figures 18 and 19 is coupled with the above discussion of total air voids analysis approach.
In Figure 20, the general snapshot of the parametric study of mean effects of strain accumulation analysis can be visualized in a compact manner.

Conclusions and Further Recommendations
Wet basis M-03 type PP fiber modification of asphalt is an efficient way to alter the mechanical properties and Marshall specimens, which are tested under repeated load creep loading by the use of UTM, have shown a considerable improvement in their service lives (5 to 12 times).In addition to this fact, a novel approach for the prediction of mechanical properties such as strain accumulation development (rutting potential) at the end of repeated load creep tests carried out by UTM utilizing GP techniques has been proposed.This approach is very important in the sense that for a certain type of asphalt mixture and for predetermined testing conditions, the strain accumulation at the end of repeated load creep tests can be estimated without carrying out destructive tests with UTM or any similar testing equipment.Therefore, the rutting potential can be explored by this means in a perfect manner as a main departure from the other published studies in this field for the very first time in the literature.The reader should bear in mind that the proposed model and parametric studies are valid for the ranges of the experimental database used for modelling.The reader must also not lose sight of that the obtained results at the end of this and similar kind of studies is valid only for similar types of aggregate sources, bitumen, aggregate gradation, mix proportioning, modification technique and laboratory conditions.Therefore it is not possible to generalise the findings that have been obtained throughout these studies for another specific case.The explicit formulation of strain accumulation on the proposed GP model is also obtained through extensive analyses and presented for further use by researchers.To obtain the main effects of each variable on strain accumulation, a wide range of parametric studies has been performed by using the GP model.As a result, the proposed GP model and formulation of the available strain accumulation development of Marshall specimens is quite accurate, fast and practical for use by other researchers studying in the experimental pavement engineering field.Wheel-tracking test, which is of course a very important instrument to determine the rutting susceptibility and potential of polymer modified bituminous mixtures can be carried out in order to determine the behaviour of PP fiber modified specimens as future research.Furthermore, by utilising GP or other soft computing methods like artificial neural networks or neuro-fuzzy techniques can be a good means for the further prediction of rutting potential of PP modified asphalt mixtures which the test results are obtained via wheel-tracking test devices and these results can also be compared and correlated with the results of the previously carried out repeated creep tests via universal testing machines.

Figure 2 .
Figure 2. Strain accumulation values for the control and M-03 type polypropylene fiber modified samples (500 ms load -500 ms rest).

Figure 3 .
Figure 3. Strain accumulation values for the control and M-03 type polypropylene fiber modified samples (500 ms load -1000 ms rest).

Figure 4 .
Figure 4. Strain accumulation values for the control and M-03 type polypropylene fiber modified samples (500 ms load -1500 ms rest).

Figure 5 .
Figure 5. Strain accumulation values for the control and M-03 type polypropylene fiber modified samples (500 ms load -2000 ms rest).

Figure 6 .
Figure 6.An example of an expression tree.

Figure 8 .
Figure 8.The genetic programming results versus actual test results for strain accumulation (rutting potential).

Figure 10 .
Figure 10.Interaction plot for rest period vs. pulse counts.

Figure 11 .
Figure 11.Interaction plot for pulse count vs. rest periods.

Figure 12 .
Figure 12.Interaction plot for specimen height vs. air voids.

Figure 16 .
Figure 16.Interaction plot for air voids (V a ) vs. rest periods.

Figure 14 .
Figure 14.Interaction plot for unit weight vs. voids filled with asphalt (V f ).

Figure 20 .
Figure 20.Whole trends for the parametric study of main effects of strain accumulation (rutting potential) analysis.

Table 1 .
Physical properties of the reference bitumen.

Table 2 .
Physical properties of coarse aggregates.

Table 3 .
Physical properties of fine aggregates.

Table 5 .
Physical properties of the 3‰ polypropylene modified bitumen samples by weight of aggregate.
a denotes polypropylene modified specimens; b denotes control specimens.

Table 7 .
Statistical details and equations of the relevant regression models.

Table 8 .
Parameters of the GEP models.

Table 9 .
Statistical parameters of testing and training sets and overall results of GP models.