**3. Developed Model Results, Performance, and Validation**

The results produced by the developed models are presented in this section, and their capability to use empirical data to estimate the dynamic modulus of asphalt mixtures is evaluated.

#### *3.1. Model Performance*

The performance is first compared with the existing predictive models; modified Witczak, Hirsch, and Alkhateeb models selected from the literature are presented in Equations (10)–(13), respectively [13,14,16]

$$\log|E^{\*}| = -0.349 + 0.754(|C\_{b}^{\*}|^{-0.0052})(6.65 - 0.032\rho\_{200} - 0.0027(\rho\_{200})^{2} + \\\\ \log|V\_{b}^{\*}|^{-2})$$

$$0.011\rho\_{4} - 0.0001(\rho\_{4})^{2} + 0.006\rho\_{3/8} - 0.00014(\rho\_{3/8})^{2} - 0.08V\_{a} - 1.06(\frac{V\_{beff}}{V\_{beff} + V\_{a}}) $$

$$+ \frac{2.558 + 0.032V\_{a} + 0.713(\frac{V\_{beff}}{V\_{brf} + V\_{a}}) + 0.0124\rho\_{3/8} - 0.0001(\rho\_{3/8})^{2} - 0.0098\rho\_{3/4}}{1 + \exp(-0.7814 - 0.5785\log|G\_{b}^{\*}| + 0.8834\log\delta\_{b})} \tag{10}$$

where |*E*∗| is dynamic modulus in psi; |*G*∗| is the binder shear modulus in psi; *δ<sup>b</sup>* is the binder phase angle in degrees; *ρ*3/<sup>4</sup> is the cumulative percent aggregate retained on the 3/4 sieve (19 mm); *ρ*3/<sup>8</sup> is the cumulative percent aggregate retained on the 3/8 sieve (9.5 mm); *ρ*<sup>4</sup> is the cumulative percent aggregate retained on the No. 4 sieve (4.75 mm); *ρ*<sup>200</sup> is the percent aggregate passing the No. 200 sieve (0.075 mm); *Va* is the percent air void in the mix; *Vbef f* is the effective asphalt content; VMA is the percent of voids in the mineral aggregate, and VFA is the percent voids filled with asphalt,

$$|E\_{\rm m}^\*| = P\_{\varepsilon} \left( 4, 200, 000(1 - \frac{VMA}{100}) + 3 \vert G^\* \vert\_{\delta} (\frac{VFA \times VMA}{10, 000}) \right) + \frac{(1 - P\_{\varepsilon})}{\frac{1 - VMA}{4, 200, 000} + \frac{VMA}{3 \vert G\_{\rm b}^\* \vert (VFA)}} \tag{11}$$

where

$$P\_{\mathcal{E}} = \frac{(20 + 3|G\_b^\*| (VFA) / (VMA))^{0.58}}{650 + (3|G^\*|\_b (VFA) / (VMA))^{0.58}} \tag{12}$$

and |*E*∗|*<sup>m</sup>* is dynamic modulus of HMA in psi; *Pc* is the aggregate contact volume; VMA is the percentage of mineral aggregate voids in compacted mixture; and VFA is the percentage of voids filled with asphalt in the compacted mixture,

$$|E\_m^\*| = 3(\frac{100 - VMA}{100})(\frac{90 + 1.45 \frac{[G^\*]\_b}{VMA})^{0.66}}{1100 + (0.13 \frac{[G\_g^\*]}{VMA})^{0.66}})|G\_g^\*|\tag{13}$$

where |*E*<sup>∗</sup> *<sup>m</sup>*|, |*G*<sup>∗</sup> *<sup>b</sup>* |, and |*G*<sup>∗</sup> *<sup>g</sup>* | (the complex shear modulus of binder in the glassy state, assumed to be 109 Pa.) are in Pa. Equation (14) shows the best reduced third-order (linear) regression model (PCR) fitting the measured response:

$$\begin{aligned} \mathcal{Y} &= c\_0 + c\_1 p c\_1 + c\_2 p c\_2 + c\_3 p c\_3 + c\_4 p c\_4 + c\_5 p c\_5 \\ &+ c\_6 p c\_1 p c\_2 + c\_7 p c\_1 p c\_3 + c\_8 p c\_1 p c\_4 + c\_9 p c\_1 p c\_5 \\ &+ c\_{10} p c\_2 p c\_3 + c\_{11} p c\_2 p c\_4 + c\_{12} p c\_2 p c\_5 + c\_{13} p c\_3 p c\_4 \\ &+ c\_{14} p c\_3 p c\_5 + c\_{15} p c\_4 p c\_5 + c\_{16} p c\_1 p c\_2 p c\_3 + c\_{17} p c\_1 p c\_2 p c\_4 \\ &+ c\_{18} p c\_1 p c\_2 p c\_5 + c\_{19} p c\_1 p c\_3 p c\_4 + c\_{20} p c\_1 p c\_3 p c\_5 + c\_{21} p c\_2 p c\_3 p c\_4 \\ &+ c\_{22} p c\_1 p c\_4 p c\_5 + c\_{23} p c\_2 p c\_4 p c\_5 + c\_{24} p c\_3 p c\_4 p c\_5 \end{aligned} \tag{14}$$

where, *c*<sup>0</sup> = 6.59; *c*<sup>1</sup> = 2.58; *c*<sup>2</sup> = 4.4; *c*<sup>3</sup> = −0.36; *c*<sup>4</sup> = 0.49; *c*<sup>5</sup> = 1.93; *c*<sup>6</sup> = −0.33; *c*<sup>7</sup> = −0.77; *c*<sup>8</sup> = −1.69; *c*<sup>9</sup> = 0.15; *c*<sup>10</sup> = −1.65; *c*<sup>11</sup> = −4.68; *c*<sup>12</sup> = 4.81; *c*<sup>13</sup> = 0.7; *c*<sup>14</sup> = −0.85; *c*<sup>15</sup> = −1.58; *c*<sup>16</sup> = −0.17; *c*<sup>17</sup> = −0.79; *c*<sup>18</sup> = 1.83; *c*<sup>19</sup> = 0.04; *c*<sup>20</sup> = 0.18; *c*<sup>21</sup> = 0.42; *c*<sup>22</sup> = 0.05; *c*<sup>23</sup> = 0.32; *c*<sup>24</sup> = 0.06. The trained three-layer ANN (PCNN) presented in Equation (9) contains the following connection weights and biases:

$$\mathbf{W}^{T} = \begin{bmatrix} -0.511 & 0.134 & 0.654 & -1.064 & -0.267 \\ -0.315 & -0.147 & -0.267 & 0.177 & -1.047 \\ -0.660 & -1.266 & 0.759 & -1.248 & -0.331 \\ -0.075 & 0.022 & 0.208 & 0.015 & 0.167 \\ -0.074 & 0.022 & 0.206 & 0.015 & 0.167 \\ 0.103 & -0.177 & 1.253 & -1.045 & 0.558 \\ 0.078 & -0.020 & -0.231 & -0.014 & -0.172 \\ 0.238 & 0.070 & -0.885 & 0.848 & 0.943 \\ 0.123 & 0.456 & -0.387 & 1.547 & -0.017 \\ -0.079 & 0.020 & 0.213 & 0.014 & 0.173 \end{bmatrix}$$

$$\mathbf{W}\_{H} = \begin{bmatrix} 0.869 \\ -0.886 \\ -0.866 \\ 0.632 \\ -0.291 \\ -0.859 \\ 0.299 \\ 0.299 \\ 0.007 \\ 0.007 \\ 0.007 \\ -0.373 \\ -0.373 \\ -0.007 \end{bmatrix}, \quad \mathbf{B}\_{0} = \begin{bmatrix} 0.162 \\ 0.150 \\ 0.250 \\ 0.070 \\ 0.070 \\ 0.071 \\ -0.373 \\ -0.007 \end{bmatrix}$$

Figure 6 presents the performance of the developed models in terms of measured values of dynamic modulus versus the fitted dynamic modulus values. The measured and fitted values are fairly close to the line of equality, indicating that the fitted values are highly correlated with the measured ones.

**Figure 6.** Measured values of dynamic modulus versus fitted values by PCR and PCNN.

Comparisons of PCR and PCNN performance to that of the existing predictive models are conducted based on three statistics: average difference (AD), average absolute difference (AAD), and correlation between measured and fitted values of response (*rfit*). A summary of the definitions of these statistical components and their formulas is presented in Table 5. In the formulas presented in Table 5, *yi* is the *i*th measured response, *y*ˆ*<sup>i</sup>* is the *i*th fitted response, and *n* is the number of data points.

The results of the comparison are presented in Table 6. According to the values of *rfit* in Table 6, the estimated dynamic modulus values obtained form PCR and PCNN models are highly correlated

with measured values according to the values, showing that the both PCR and PCNN performed well in terms of modeling the response variable.


**Table 5.** Statistics which are used to compare model performance.

Although the corresponding values of *rfit* for modified Witczak, Hirsch, and Alkhateeb models are 0.93, 0.95, and 0.95, respectively, the average difference and average absolute difference with respect to the measured response are significantly higher than those of PCR and PCNN. This means that the fitted values by the modified Witczak, Hirsch, and Alkhateeb models are not close as those fitted by PCR and PCNN to the response value. In other words, *rfit*, which reflects the correlation between response and estimated response (if one goes up the other one goes up), could be biased, and in this situations other statistics (AD, and AAD) could be used to evaluate the goodness of fit. The dynamic modulus measured and predicted values are presented in Figure 7 for four asphalt mixtures. According to the presented master curves the current study (PCNN model) provides the closest values of *E*∗ to the measurements for all of three test temperatures, while, the conventional models either overestimate or underestimate the response variable.

**Figure 7.** Comparing the measured and predicted (current study and conventional models) dynamic modulus. Results are presented at three temperatures. Green color indicates 0.4 ◦C, Orange color indicates 17.1 ◦C, and Purple color indicates 33.8 ◦C.



A graphical comparison of the PCR and PCNN performance and that of the existing models is presented in the following section.

#### *3.2. Receiver Operating Characteristic Analysis (ROC)*

A receiver operating characteristic (ROC) graph is a technique for visualizing, organizing, and selecting classifiers based on their performance. ROC graphs are widely used in medical decision-making as well as in machine learning and data-mining research [35]. True ROC curves plot the false positive rate (probability of false alarm) on the *x*-axis and the true positive rate (probability of detection) on the *y*-axis. A classifier is said to perform well if the ROC curve climbs rapidly towards the upper left-hand corner. The more the curve deviates from *y* = *x* behavior, the more accurate the prediction is [36]. We can borrow from the concept of ROC curve to obtain a measurement of fit for the competing models and a ROC graph for this study is presented in Figure 8 for this study. As described in Equation (15), the *x*-axis indicates the standardized residuals ordered from the lowest to the highest (*e*∗∗ *<sup>i</sup>* ). Residuals are sorted in ascending order and divided by the largest one that belongs to the Hirsch model.

$$\left|e\_{i}^{\ast \ast}\right| = \frac{\left|e\_{i}^{\ast}\right|}{\left|e\_{i\_{M\mathcal{A}\mathcal{A}}}^{\ast}\right|}\tag{15}$$

**Figure 8.** ROC curves for the developed and the existing predictive models for the dynamic modulus. The curves indicate the goodness of the fit provided by the models with regard to a pre-specified residual. The PCNN model showed the highest performance (farthest from the *y* = *x* line), and the Hirsch model showed the lowest (closest to the *y* = *x* line).

The *y*-axis indicates the fraction of points whose standardized residuals are less than *e*∗∗ *<sup>i</sup>* [37]. Although the curves obtained for all of the models are monotonically non-decreasing and climb towards the upper left-hand corner (the desired situation shows that the predictive models perform well), PCNN and PCR curves were the highest, proving their better performance to be better than that of existing predictive models.

A convenient global measure of the goodness-of-the-fit is the Area Under the Curve (AUC). To compare classifiers, it is more desirable to reduce ROC performance, to a single scalar value representing expected performance. Since the AUC is a portion of the area of the unit square, its value will always lie between 0 and 1. The AUC values for the PCNN, PCR, Alkhateeb, modified Witczak, and Hirsch models are 0.9864, 0.9717, 0.8746, 0.7609, and 0.6320, respectively. One can use the ROC and AUC analysis results and rank the predictive models according to their performances. In this study, the PCNN model reflected the highest performance in predicting the dynamic modulus value, while the Hirsch model ranked the lowest among all the models.

#### *3.3. Model Validation*

The current regression model is presented in the following general form as

$$y\_i = f\_i(\mathbf{Z}\_i, \boldsymbol{\theta}) + \boldsymbol{\varepsilon}\_i^\*. \tag{16}$$

In the above equation *fi* is the *i*th expectation function, *θ* is the vector of parameters, and *e*<sup>∗</sup> *<sup>i</sup>* is a random deviation of *yi* from *fi*. This term is assumed to be independent and normally distributed with a mean of zero and unknown variance *<sup>σ</sup>*<sup>2</sup> for *<sup>i</sup>* = 1, ···, *<sup>n</sup>*, where *<sup>n</sup>* is the number of input vectors. If the above assumptions are violated, the results of the analysis could be misleading or erroneous. These assumptions can be testified by examining residuals as defined by

$$e\_i^\* = \mathcal{Y}\_i - \mathcal{Y}\_i. \tag{17}$$

The assumption of independency holds when the residuals plot does not reflect a trivial pattern. The normality assumption is assessed by creating a normal probability plot of the residuals. When the error has a normal distribution, this plot will appear as a straight line [10]. These assumptions were checked for PCR and PCNN, as presented in Figure 9. The assumption of equal variances does not appear to be violated because there are no trivial pattern in this plot. Figure 9 presents the normal probability of the residuals in which it can be seen that the data points are close to the straight line and the normality assumption is validated.

**Figure 9.** *Cont.*

**Figure 9.** Checking the assumptions of independency using residual plot for PCR (**left**) and PCNN (**right**) models and the normality using normal probability plot of the residuals for PCR (**left**) and PCNN (**right**) models.

#### **4. Application of the Framework: Flexible Pavement Design and Optimization**

The above framework is used along with an optimization algorithm to answer the following two central questions:


One can see that the first item corresponds to the optimal design problem while the second one corresponds to the so-called inverse design.

Since it was shown through multiple statistical measurements that PCNN had the best prediction capability, this model is used in the following section to solve the optimization problems. The ANN used in PCNN is essentially an interconnected nonlinear function, and this necessitates the application of a global optimizer. Moreover, the effective variable space enters the problem as a series of constraints and further restricts the available algorithms. The optimal design problem is formulated as follows:

$$\begin{aligned} \text{maximize} \quad & |E^\*| = F\_{ANN}(\mathbf{x})\\ \text{with respect to} \quad & \mathbf{x} = (\mathbf{x}\_1, \dots, \mathbf{x}\_{14})\\ \text{subject to} \quad & (\mathbf{x} - \mathbf{v})^T A (\mathbf{x} - \mathbf{v}) \le \mathbf{1},\\ & (\mathbf{x}\_{pca} - \mathbf{v}')^T A' (\mathbf{x}\_{pca} - \mathbf{v}') \le \mathbf{1}, \end{aligned} \tag{18}$$

where the vector of fourteen variables is **x**, and (**x** − **v**) *TA*(**<sup>x</sup>** − **<sup>v</sup>**) ≤ **<sup>1</sup>** are the enclosing ellipsoid constraint equations for the original and PCA-based variables. A penalty function approach is used to convert the above constrained problem to an unconstrained one [38]. In this case, when the penalty function is active, it decreases (increases) the objective function when the problem is one of maximization (minimization), and the degree of penalty is based on the closeness of the solution to the corresponding constraint.

Since the inverse design problem aims at finding the specification of a predefined goal, it is defined as a minimization problem as follows:

$$\begin{aligned} \text{minimize} \quad &error = ||E^\*|| - |E\_0^\*|| \\ \text{with respect to} \quad & \mathbf{x} = (\mathbf{x}\_1, \dots, \mathbf{x}\_{14}) \\ \text{subject to} \quad & (\mathbf{x} - \mathbf{v})^T A (\mathbf{x} - \mathbf{v}) \le \mathbf{1}, \\ & (\mathbf{x}\_{pca} - \mathbf{v}')^T A' (\mathbf{x}\_{pca} - \mathbf{v}') \le \mathbf{1}, \end{aligned} \tag{19}$$

where |*E*<sup>∗</sup> <sup>0</sup> | is the desired (goal) dynamic modulus. Although a similar penalization method can also be used to address the constraints in this case, for the above problem the constraints will penalize the objective function when they are active.

Reliable solution of the above problems requires the application of a gradient-free optimization algorithm. Gradient-based optimization algorithms are not applicable in this case because of the network-based nature of the ANNs. Evolutionary-based algorithms are potentially easy-to-use algorithms in the above problems. Novel algorithms have been used to solve complex optimization problems in recent years [39,40], and in this case, Mean-Variance Mapping Optimization (MVMO), an in-house optimization algorithm based on the work by Elrich et al. [41,42], is used. The constraints are handled using the approach described in Aslani et al. [43], in which, the convergence rate of a constrained MVMO was compared to the already-developed methodologies using benchmark structural problems. Authors in [22] indicated that a constrained MVMO is capable accurately identifying an optimal value with a minimum number of simulations. It should be noted that the choice of optimization algorithm is not the principal focus of this study.

Figure 10 (left) depicts the convergence achieved for the first design problem by the constrained MVMO algorithm. The initial data points are random making it heavily penalized, and then the objective function increases as the algorithms evolves. Exploration-exploitation behavior is achieved using adaptive strategies in the course of optimization for MVMO. *δ* = 0.05 is used as the threshold in Figure 4. Solving the maximization problem resulted in |*E*<sup>∗</sup> *max*|= 53,703 MPa. The optimal design parameters are presented in the first column of Table 7.

To find the maximum amount of dynamic modulus one could design for without low temperature failure in the asphalt binder, the maximization problem was solved one more time with an additional constraint of *G*∗*sinδ* ≤ 5000, resulting in |*E*<sup>∗</sup> *max*|= 36,307 MPa. Corresponding design parameters are presented in the second column of Table 7 as the optimal design 2.

Figure 10 (right) shows the convergence of the algorithm for the inverse design problem after starting randomly from three different initial points, with the algorithm is terminated when the error reaches about 10<sup>−</sup>9. A pre-specified |*E*<sup>∗</sup> <sup>0</sup> | of 20,417 MPa is considered and the inverse problem of finding the corresponding design parameters is solved. Because of non-linearity of the function, the problem has no unique solution. Three of the possible solutions are presented as designs 1 to 3 in Table 7.

Finally, the five sets of design parameters are compared with current design specification, with the results shown in Table 7. The percentage of aggregate passing by each sieve size is within the acceptable range of the gradation specification. Gradation charts are presented in Figure 11. The obtained percentages of air voids are 4%, which is the target value in the design specification. The obtained values for VMA are slightly less than 14% for a nominal maximum aggregate size (NMAS) of 12.5 mm because the VMA values of the nine mixtures used to train the PCNN are slightly less that 14% (see Table 1). The acceptable range for VFA varies with the amount of traffic load measured in million Equivalent Single Axle Loads (ESALs) as follows:


The VFAs obtained for all of the five sets of design are satisfied for all of the traffic categories.

**Figure 10.** Convergence for the optimal design problem (**left**) and inverse (**right**) design of dynamic modulus.

**Table 7.** Th Corresponding Design Parameters Obtained from Solving Optimal Design and Inverse Design Problems.


**Figure 11.** Aggregate gradation graphs with 12.5 mm NMAS particle size distribution obtained from PCNN.

#### **5. Conclusions**

This study used the HMA dynamic modulus data and focused to evaluate the quality of predictor variables to be used in a procedure of model development. Correlation analysis is performed to identify cross-correlated input variables, and correlated inputs are replaced by orthogonal pseudo-inputs (PCs) obtained using PCA. Two separate models are developed using multivariate regression and ANN (called PCR and PCNN, respectively). Extrapolation in empirical modeling is addressed by adding the constraint of an n-dimensional enclosing ellipsoid to the modeling problem. Performances of the proposed models were compared to existing predictive models using both statistical analysis and ROC analysis. The models developed satisfactorily estimated the dynamic modulus value, with PCNN indicating remarkably better performance when fitted to the test data than the existing predictive models from the literature. These PCA-based approaches are thus highly recommended as precise modeling strategies in this application. Moreover, these methodologies appear to be capable of modeling other material properties and future investigation in this regard is recommended. To determine this framework's application in pavement design, two optimization problems including optimal design and inverse design have been presented and solved using a mean-variance mapping optimization algorithm. The results for the two problems are in a good agreement with the HMA mix design specification and thus could be a reasonable starting point in solving real-life design problems. Although, the developed models as well as obtained optimal design parameters are based on the empirical database created in this study, the suggested framework has the capability of being re-trained and adjusted to fit new data. For obtaining more reliable and applicable results, a larger empirical database would be required.

**Author Contributions:** Conceptualization, P.G. and M.A.; methodology, P.G., M.A., and D.K.R.; software, P.G. and M.A.; validation, P.G., M.A., D.K.R., and R.C.W. formal analysis, P.G., M.A., D.K.R., and R.C.W.; investigation, P.G., M.A., D.K.R., and R.C.W.; resources, P.G., M.A., D.K.R., and R.C.W.; data curation, P.G., M.A., D.K.R., and R.C.W.; writing–original draft preparation, P.G., M.A., D.K.R., and R.C.W.; writing–review and editing, P.G., M.A., D.K.R., and R.C.W.; visualization, P.G., M.A., D.K.R., and R.C.W.; supervision, P.G., M.A., D.K.R., and R.C.W.; project administration, R.D.W.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to thank anonymous reviewers for their fruitful comments on the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.
