Parameter Determination of the 2S2P1D Model and Havriliak–Negami Model Based on the Genetic Algorithm and Levenberg–Marquardt Optimization Algorithm

This study utilizes the genetic algorithm (GA) and Levenberg–Marquardt (L–M) algorithm to optimize the parameter acquisition process for two commonly used viscoelastic models: 2S2P1D and Havriliak–Negami (H–N). The effects of the various combinations of the optimization algorithms on the accuracy of the parameter acquisition in these two constitutive equations are investigated. Furthermore, the applicability of the GA among different viscoelastic constitutive models is analyzed and summarized. The results indicate that the GA can ensure a correlation coefficient of 0.99 between the fitting result and the experimental data of the 2S2P1D model parameters, and it is further proved that the fitting accuracy can be achieved through the secondary optimization via the L–M algorithm. Since the H–N model involves fractional power functions, high-precision fitting by directly fitting the parameters to experimental data is challenging. This study proposes an improved semi-analytical method that first fits the Cole–Cole curve of the H–N model, followed by optimizing the parameters of the H–N model using the GA. The correlation coefficient of the fitting result can be improved to over 0.98. This study also reveals a close relationship between the optimization of the H–N model and the discreteness and overlap of experimental data, which may be attributed to the inclusion of fractional power functions in the H–N model.


Introduction
Polymers are widely used in various fields, such as construction materials, electronic devices, high-energy explosives, textiles, and medical applications due to the fact of their excellent mechanical properties, including high deformability, damping characteristics, and strong adhesion to other materials [1][2][3][4]. Since World War II, polymeric materials have come to occupy a position of equal importance to metallic materials in the modern industrial system.
In general, polymers exhibit significant temperature-and time-dependent viscoelastic behavior, which is determined by the molecular structure of polymer compounds. Unlike metallic materials, polymers are formed by the polymerization of many small molecules into large polymers, resulting in an interconnected network structure. The viscoelasticity of polymer materials is caused by the deformation of the chain-like polymers under stress. With the widespread use of polymer materials in various engineering fields, an accurate assessment of the viscoelastic properties of polymers has become crucial for predicting their mechanical behavior [5,6].
To accurately describe the viscoelastic response of polymer materials, various constitutive models have been proposed [7]. Common viscoelastic models include the Maxwell model, Kelvin-Voigt model, generalized Maxwell model, standard linear model (SLM), 2S2P1D model, and Huet model [8]. Among these, the 2S2P1D model is a generalized model derived from the Huet-Sayegh model that accurately describes the rheological properties of adhesives and asphalt mixtures [9]. The Havriliak-Negami (H-N) model is of great value in predicting the viscoelastic properties of rubber-like materials and their engineering applications due to the fact of its simplicity and accurate prediction capabilities [10].
Different viscoelastic constitutive models require the determination of different types of viscoelastic parameters. Since viscoelastic constitutive models involve time information, it is difficult to express their expressions directly. They are usually represented using a discretization method. Common discretization methods for viscoelastic constitutive model parameters include Prony series models and Schapery constitutive models [11,12]. Since the parameters of viscoelastic constitutive models are difficult to measure directly, they need to be determined by comparing optimization methods with simulation results [13][14][15][16]. Since Prony series and Schapery constitutive models have a large research base and open source programs for solving the time-varying behavior of polymer materials and structures in complex states in combination with finite element methods, they are widely used for solving the time-dependent behavior of complex states in polymer materials and structures at present [17][18][19][20]. However, the description of the viscoelastic behavior of materials is usually achieved by using various functions, such as exponential, sigmoid functions, and various Bessel functions, to define continuously varying elastic moduli for different loading periods, which are mathematically represented in the form of an integral equation with a graded integral kernel. Mathematically, this is represented by an integral equation containing a series integral kernel [21][22][23][24]. In practical applications, the integral equation needs to be discretized, which means discretizing the continuous exponential function into Prony series. The number of sampling points in the discretization process determines the computational accuracy, and the number of sampling points in the Prony series determines the number of parameters in the constitutive model. Having too many sampling points obviously increases the difficulty of determining the parameters and may lead to overfitting issues [25][26][27][28]. To solve this problem, extensive research has been conducted over the past two decades on methods related to the identification and acquisition of viscoelastic parameters using various optimization algorithms, including gradient-based batch gradient descent, stochastic gradient descent, mini-batch gradient descent, and the Levenberg-Marquardt (L-M) algorithm, as well as gradient-free optimization algorithms, such as genetic algorithm (GA), particle swarm optimization, and neural network optimization algorithms. There are two common methods for obtaining viscoelastic parameters in polymer materials: one is the fitting method, based on relaxation and creep test results, and the other is the fitting method, based on the master curve of the dynamic modulus. The method based on dynamic modulus testing is favored by researchers due to the fact of its high accuracy and wide frequency range. Although gradient and non-gradient optimization algorithms have been widely used in the inverse analysis of viscoelastic model parameters, there have been few mixed optimization algorithms that combine these two methods for the parameter acquisition and optimization of different materials and viscoelastic constitutive models. The use of the GA combined with the nongradient optimization algorithm to obtain the initial values of the parameters, as well as utilizing the L-M algorithm, a gradient algorithm, for the further optimization of the parameters, is of great significance for improving the accuracy of viscoelastic constitutive model parameters and studying the influence of viscoelastic materials in complex mechanical time-dependent processes.
Based on this, this paper intended to explore the optimization problem of viscoelastic model parameters using different optimization algorithms. Section 2 of this paper mainly introduces the basic theory of GA and their applications in the parameter acquisition and optimization of the 2S2P1D model and H-N model. Detailed calculation processes for inverting parameters are provided. In Section 3, based on experimental data, two methods of optimization inversion are studied for the 2S2P1D model and H-N model: one is based on empirical parameters, and the other is based on the William, Landel, and Ferry (WLF) equation of relaxation time. Furthermore, based on the preliminary fitting results using the GA, a sequential mixed algorithm using GA+L-M is proposed to further optimize the parameters of the two models, and the influence of the L-M algorithm on the optimization results of the GA is investigated. Finally, the impact of the GA and GA+L-M algorithms on obtaining parameter results in these two commonly used viscoelastic constitutive equations is discussed, and the applicability of the GA and GA+L-M algorithms among different models is summarized.

Methodology and Optimization Study
In this section, two commonly used viscoelastic models, including the 2S2P1D and H-N models, are applied. The formula of the H-N model contains high-order fractional derivatives, complex numbers, and generalized integrals, which place high demands on the selected upper and lower bounds of parameters. Thus, the direct optimization of the dynamic modulus and phase angle usually results in inaccurate values or even a large discrepancy between the solution and the fitting results due to the fact of inappropriate initial values. Figure 1 illustrates the basic principles of the GA. It is a method that searches for the global optimal solution by simulating the natural evolution process. The algorithm begins with an initial population and utilizes random selection, crossover, and mutation operations to generate individuals that are better adapted to the environment. This iterative process promotes continuous evolution, leading to convergence on a group of individuals highly adapted to the environment, ultimately obtaining the optimal solution to the problem at hand [29]. The fundamental components and operations of the GA include the following: (1) Population-a GA maintains a large population of individuals, representing candidate solutions to the problem. Each individual is represented by a chromosome, and the population can be viewed as a collection of chromosomes. (2) Genotype-in the GA, each individual consists of a chromosome representing a set of genes. (3) Fitness function-at each iteration of the algorithm, individuals are evaluated using a fitness function (also known as an objective function), which quantifies the quality of a solution. (4) Selection-After calculating the fitness of each individual, a selection process determines which individuals should be used to breed and produce the next generation. (5) Crossoverto create new individuals, portions of parent chromosomes from the current generation are exchanged, resulting in two descendant chromosomes. (6) Mutation-The mutation operation simulates gene mutations occurring in the natural genetic environment. It randomly changes the value of a gene with a small probability. In binary encoding, mutation manifests as the random flipping of a gene from 1 to 0 or vice versa. The pseudocode for the GA implemented in this study is provided in Appendix A.

Theory of 2S2P1D Model
The 2S2P1D model, initially calibrated at the Ecole Nationale des Travaux Publics de l'Etat (ENTPE) laboratory [30], is a widely used viscoelastic constitutive model for modeling the behavior of viscoelastic materials. "2S2P1D" stands for the abbreviation of two springs, two parabolic creep elements, and one damper. As shown in Figure 2a

Theory of 2S2P1D Model
The 2S2P1D model, initially calibrated at the Ecole Nationale des Travaux Publics de l'Etat (ENTPE) laboratory [30], is a widely used viscoelastic constitutive model for modeling the behavior of viscoelastic materials. "2S2P1D" stands for the abbreviation of two springs, two parabolic creep elements, and one damper. As shown in Figure 2a, the model is an extension of the Huet-Sayegh model. Its basic theory is as follows:  Based on research of Olard and Di Benedetto, as shown in Reference [9], the parabolic creep element obeys the parabolic creep function: The expression of the complex modulus is: where J(t) is the creep function; h is the exponential index; 0 < h < 1(h = 0 (elastic), h = 1(viscous)); a is a dimensionless constant; λ is the relaxation time; Γ is the gamma function; t is the loading time; τ is the relaxation time (varies only with temperature); i is an imaginary number; and ω is the frequency. For the linear damper, the complex modulus can be expressed as follows: The expression of the 2S2P1D model is shown as follows: where k and h are the exponent indices, with 0 < k < h < 1; α is a constant; G 0 is the static modulus at ω → 0 ; and G g is the glassy modulus at ω → 1 . As shown in Figure 2b, G 0 and G g are the intercepts of the Cole-Cole curve on the x-axis. β is a constant. η is defined as: where η is the Newtonian viscosity, and τ is the relaxation time, which is a function of the temperature, as shown in Figure 2c. Within the temperature range observed in the laboratory, τ can be approximated by the shift factor law based on the WLF equation: where C 1 and C 2 are empirical parameters, T and T r are the actual and reference temperatures, a T is the shift factor, and τ 0 and τ i are the relaxation times at temperatures T r and T, respectively. The WLF equation is a significant empirical formula in polymer physics that delineates the correlation between the relaxation time and temperature. In numerous amorphous polymers, the complex modulus at different temperatures can be superimposed into a master curve with a broader frequency range by horizontally shifting the time quantities obtained at different temperatures. The phase angle can be expressed as: where G and G are the storage and loss moduli. From Equation (4), we can see that seven parameters (G 0 , G g , k, h, α, β, and τ) are required to be determined in the 2S2P1D model.

Parameter Optimization Process
In this section, the GA and L-M algorithms are used to conduct the optimization of the 2S2P1D model's parameters. Two methods (dependent and independent WLF functions) for the optimization are considered. The optimization processes are as follows: (1) Parameter optimization independent of the WLF function In the first approach, the test data obtained at different temperatures are independently optimized, resulting in the determination of the relaxation time (τ i ) values for each temperature. Based on Equation (4), we know that six temperature-independent parameters (G 0 , G g , k, h, α, and β) and one temperature-dependent parameter (τ i ) should be determined. The detailed procedure is given below: (1) Take G 0 , G g , k, h, α, β, and τ i as unknowns and the regularized 2-norm sum of the differences between the calculated (based on Equations (4) and (8)) and experimental dynamic moduli and phase angles as the objective function, and then conduct the optimization using the GA considering the given constraints and upper and lower bounds. (2) Take the parameters obtained from the GA as the initial values for further optimization using the L-M algorithm.
Substitute the obtained parameters into Equations (4) and (8) to calculate the corresponding dynamic modulus (|G * |) and phase angle (δ). The correlation coefficient between the calculated values and the original experimental values is calculated to evaluate the fitting results.
(2) Parameter optimization dependent on WLF function In the second method, the test data obtained at different temperatures are shifted to the reference temperature (Tr) using the WLF equation, with the subsequent optimization of only one relaxation time value (τ 0 ) at the reference temperature. The specific steps are as follows: (1) Take G 0 , G g , k, h, α, β, τ 0 , C 1 , C 2 and T r as unknowns, and employ the GA to obtain the initial parameter values. (2) Take the parameter values obtained with the GA as the initial values and further optimize the parameters using the L-M algorithm.
Substitute the obtained parameters into Equations (4) and (8) to calculate the corresponding dynamic modulus (|G * |) and phase angle (δ). The correlation coefficient between the calculated values and the original experimental values is calculated to evaluate the fitting results.

Theory of the H-N Model
The H-N model is widely used to characterize the complex dielectric behavior of polymers. The relationship between the complex modulus and the instantaneous modulus (E 0 ) and the long-term modulus (E ∞ ) given by the equation is [10]: where ω is the angular frequency; τ is the relaxation time related to temperature; i is the unit imaginary number; E 0 , E ∞ , α, and β are the parameters of the H-N model that are independent of temperature; and E (ω) and E (ω) are the storage modulus and loss modulus, respectively. However, the H-N model is highly intricate and involves fractional and nonreal power functions. This complexity can lead to nonreal outcomes if the parameter analysis is directly performed on these functions. To address this concern, this paper employed a semi-analytical method to conduct the inverse analysis on the constitutive model. This approach involves plotting the Cole-Cole curve and utilizing the semianalytical method to determine the polynomial coefficients. Subsequently, these coefficients are substituted into the Cole-Cole diagram using the four tangent equations, which are presented below: The temperature-independent parameters of the H-N model, namely, the four angles (ϕ, θ, Φ, and Θ), can be initially determined based on Equations (10)-(15). In Figure 3, the Cole-Cole diagram is illustrated. Typically, a Cole-Cole diagram exhibits an inverted U-shaped curve for viscoelastic materials. The OE axis is the storage modulus axis, and the OE axis is the loss modulus axis. The curve E 0 AE ∞ is the Cole-Cole curve, which is typically employed to depict the relationship between the loss modulus and storage modulus of viscoelastic materials [10].
Ê , αˆ, and βˆ, determine the optimization interval and apply the genetic algorithm again to optimize the temperature-independent parameters ( 0 E , ∞ E , α , and β ) based on the objective function, as shown below: (20) where exp E′ and exp E′′ represent the experimental storage modulus and loss modulus, respectively. The expressions of E′ and E ′ ′ are presented in Appendix C for Equations (A3) and (A12), respectively.

Parameter Optimization Process
Previous studies have directly optimized the H-N model's parameters based on the dynamic modulus and phase angle, a method highly sensitive to the optimization range and initial values [10]. Moreover, the H-N model formula incorporates high-order fractional derivatives, complex numbers, and generalized integrals, which impose strict requirements on the selection of the upper and lower parameter bounds. Consequently, the direct optimization of the dynamic modulus and phase angle often yields inaccurate values and may result in substantial discrepancies between the solution and the fitting results due to the fact of inappropriate initial values. In order to solve these issues and ensure numerical stability in the calculation of the complex H-N model, we propose a progressive method.
In contrast to the direct optimization approach, the Cole-Cole curve, which describes the relationship between E and E , is initially fitted using a third-order polynomial function to determine certain parameter values of the H-N model. Subsequently, the GA is employed to further optimize the parameter values. The detailed optimization process based on the proposed progressive method is outlined below: (1) Plot the storage modulus and loss modulus obtained from experimental tests on the same coordinate system. The inverted U-shaped Cole-Cole curve formed by the scattered points is approximated using a third-order polynomial, as follows: where p 1 , p 2 , p 3 , and p 4 are the coefficients of the third-order polynomial. Their values can be determined by fitting this polynomial to the test points based on the principle of least mean squares (L-M-S). (2) Based on the determined Cole-Cole curve using the polynomial, calculateÊ 0 ,Ê ∞ , dE /dE , tan ϕ, and tan θ. (3) Calculateα andβ using Equation (10) and determine the peak point of the Cole-Cole curve (E (ω n ), E (ω n )). Substitute this point into Equation (17) to check if (ω n τ)α is smaller than 1: (4) Use the GA to calculateα andβ. If (ω n τ)α < 1, the objective function can be written as Equation (18): If (ω n τ)α > 1, the objective function can be written as Equation (19).
(5) Based on the calculatedÊ 0 ,Ê ∞ ,α, andβ, determine the optimization interval and apply the GA again to optimize the temperature-independent parameters (E 0 , E ∞ , α, and β) based on the objective function, as shown below: where E exp and E exp represent the experimental storage modulus and loss modulus, respectively. The expressions of E and E are presented in Appendix C for Equations (A3) and (A12), respectively. (6) Based on the optimized values of the four parameters (E 0 , E ∞ , α, and β), apply the GA again to calculate the relaxation time (τ i ) at each temperature. Based on the WLF equation, calculate the shift factor α T (T) at each test temperature with respect to the reference temperature (Tr = 30°C). α T (T) defines the relationship between τ i and τ re f , as shown below: According to the shift factors α T (T) at different temperatures, the reduced frequency can be calculated based on Equation (22), as shown below: By shifting the experimental data at different temperatures, construct the master curve with a wider frequency range at the reference temperature. Additionally, input the parameters (E 0 , E ∞ , α, β, and τ i ) obtained through inverse analysis into the H-N model to calculate the E and E at different temperatures and plot their master curves. Finally, compare the experimental data-based master curves with those constructed using the optimized parameters.

Case Study on the Optimization of the 2S2P1D Model
In this case study, dynamic test data of asphalt were utilized for optimizing the parameters of the 2S2P1D model. The dynamic shear rheometer (DSR) was employed to conduct tests at five temperatures (ranging from 0 to 40°C) and frequencies from 0.1 Hz to 100 Hz.

Parameter Optimization Independent of the WLF Function
To optimize the parameters, the experimental data of the dynamic modulus |G * | were read, and G 0 , G g , k, h, α, β and five different relaxation times (τ i ) corresponding to five temperatures were used as the unknowns. Then, the model parameters were optimized using the GA. According to the parameters' physical meaning, the following constraints should be satisfied: These constraints facilitate the optimization process of the GA. The optimized values were obtained by minimizing the objective function defined in Equation (24), where |G * | exp and δ exp represent the experimental dynamic modulus and phase angle, respectively. The expressions for the dynamic modulus (|G * |) and phase angle (δ) are given by Equations (4) and (8), respectively.
The parameter values for the GA solver are shown in Table 1. The fitted values of the parameters for the 2S2P1D model are presented in Table 2. Furthermore, the values were further optimized using the L-M algorithm, as shown in Table 3.       Using the results obtained by applying the GA as the initial values, the L-M algorithm was employed for further optimization. By substituting the parameter values obtained from the first and second optimizations into Equations (4) and (8), we calculated the fitted dynamic modulus and phase angles. Additionally, we calculated the correlation coefficients between the experimental data and the values obtained from the GA. Figure 4a,b,d,e illustrate the experimental data, the calculated data obtained from the GA, and the quadratic fitted on the same coordinate system. Figure 4c,f depict the correlation coefficients. The pseudocode for the implementation of the L-M algorithm in this paper can be found in Appendix B and the parameter values for the L-M algorithm solver are shown in Table A1.
The optimized values by inverse calculation based on the GA were taken as the initial values and input into the L-M algorithm for secondary fitting. The general pseudocode for implementing the L-M algorithm in this paper is provided in Appendix B. The optimized values were then used in the 2S2P1D model (Equation (4)) to obtain the calculated dynamic modulus and phase angle. Figure 4a,b,d,e present the experimental data and the corresponding calculated curves. Figure 4c,f display the correlation coefficients. The pseudocode for implementing the L-M algorithm in this paper is included in Appendix B.

Parameter Optimization Dependent on the WLF Function
The experimental data of the dynamic modulus (|G * |) and phase angle (δ) were used as inputs, and the parameters G 0 , G g , k, h, α, β, τ i , C 1 , C 2 and T r were considered for optimization. These parameters have specific physical interpretations, and the following relationships should hold: The constraints were applied, and the upper and lower bounds were set for the GA. The optimized parameter values were determined by minimizing the objective function defined in Equation (26). To avoid overfitting, the data were normalized. The GA was configured with the parameter values shown in Table 4. The optimization process took 136.74 s. The resulting optimized parameter values for the 2S2P1D model are presented in Table 5.     angle are more sensitive to the selection of the initial values when fitting the physical quantities of the viscoelastic properties of asphalt.  The secondary optimization is performed using the L-M algorithm, and the resulting fitted parameter values are shown in Table 6. The calculation process took 2.15 s. By substituting the optimized parameters obtained from the L-M algorithm into the 2S2P1D model, the dynamic modulus and phase angle are calculated. Figure 4g,h,j,k illustrate the experimental data and the corresponding calculated curves, while Figure 4i,l depict the correlation coefficients. Table 6. Optimized parameters' values based on the L-M algorithm. As depicted in Figure 4c,i, the correlation coefficient of the dynamic modulus shows that the accuracy of the parameters obtained from both the GA and the L-M algorithm improved with increasing temperature. This suggests that the 2S2P1D model is more suitable for describing the dynamic modulus of asphalt at higher temperatures. Moreover, the correlation coefficients between the calculated values from the fitted models and the experimental data were consistently above 0.998, regardless of whether single fitting or secondary fitting was performed. It is noteworthy that Figure 4i demonstrates an improvement in the secondary fitting compared to the results obtained solely from the GA as the temperature decreased. This highlights the accuracy and practicality of the GA in finding optimal solutions.
In Figure 4f,l, the correlation coefficient of the phase angle reveals an enhancement in the accuracy of the fitting results obtained from both the GA and the L-M algorithm as the temperature decreased. This indicates that the 2S2P1D model provides more precise descriptions of the phase angle of asphalt under lower temperature conditions. The correlation coefficients between the results obtained from both the single fitting and secondary fitting and the experimental data are above 0.820. However, as the temperature decreased, the correlation coefficient of the secondary fitting gradually approached 1. In contrast to the dynamic modulus, the utilization of the second-order fitting consistently yielded better results in the phase angle correlation coefficient histogram compared to the use of the GA alone. This implies that the fitting results of the phase angle are more sensitive to the selection of the initial values when fitting the physical quantities of the viscoelastic properties of asphalt.
In this study, the parameters obtained from the GA were incorporated into the 2S2P1D model. Numerical values within a predetermined range were generated and input into Equations (4) and (8) of the 2S2P1D model to calculate the corresponding dynamic modulus and phase angles. The calculated curves were then compared with the experimental data curves. Additionally, correlation coefficients were calculated between the experimental data and the values calculated using the GA and the secondary fitting using the L-M algorithm, as depicted in Figure 5. The results show that the trend of the modeled curves was consistent with the experimental data, indicating satisfactory fitting of the dynamic modulus and phase angle. The optimization methods successfully determined suitable values for the 2S2P1D model, which effectively describes the viscoelastic properties of asphalt.

Case Study on the Optimization of the H-N Model
The test data used in this case study were derived from the dynamic mechanical thermal analysis (DMTA) data of carbon-black-filled rubber (CBFR) NR100N550 [31]. By fitting the experimental data points of the Cole-Cole curve with a third-order polynomial function, the parameters (α, β, and γ) of the polynomial function (see Equation (17)) could be determined, as shown in Table 7.
namic modulus and phase angles. The calculated curves were then compared with the experimental data curves. Additionally, correlation coefficients were calculated between the experimental data and the values calculated using the genetic algorithm and the secondary fitting using the L-M algorithm, as depicted in Figure 5. The results show that the trend of the modeled curves was consistent with the experimental data, indicating satisfactory fitting of the dynamic modulus and phase angle. The optimization methods successfully determined suitable values for the 2S2P1D model, which effectively describes the viscoelastic properties of asphalt.

Case Study on the Optimization of the H-N Model
The test data used in this case study were derived from the dynamic mechanical thermal analysis (DMTA) data of carbon-black-filled rubber (CBFR) NR100N550 [31]. By fitting the experimental data points of the Cole-Cole curve with a third-order polynomial function, the parameters (α, β, and γ) of the polynomial function (see Equation (17)) could be determined, as shown in Table 7. Table 7. Fitted parameters of the third-order polynomial function. By substituting the experimental data of into the polynomial, the corresponding " could be obtained [10]. Based on Equation (17), setting = 0, the corresponding two estimated values, which correspond to 0 E and ∞ Ê , could be obtained. With further   By substituting the experimental data of E into the polynomial, the corresponding E could be obtained [10]. Based on Equation (17), setting E = 0, the corresponding two estimated E values, which correspond toÊ 0 andÊ ∞ , could be obtained. With further calculations of the derivative values of dE /dE atÊ 0 andÊ ∞ , the corresponding values of tan ϕ and tan θ could be derived. The results are presented in Table 8. Table 8. EstimatedÊ 0 ,Ê ∞ tan ϕ and tan θ. According to Equation (11), the estimations were as follows:α = 0.564 andβ = 0.400. Based on the peak point of the Cole-Cole curve (E (ω n ), E (ω n )) = (1.746 × 10 3 , 9.370 × 10 2 ). Substituteα andβ into Equation (17) to check whether (ω n τ)α < 1. Using Equation (12), the estimated values ofα andβ were calculated by optimizing the objective function with the GA, as shown in Equation (18). The lower limit of the optimization interval was [α β ](1 − ∆ 1 ), and the upper limit was [α β ](1 + ∆ 1 ). The calculation results areα = 0.606 andβ = 0.502. With the obtainedÊ 0 ,Ê ∞ ,α, andβ, the optimization interval was determined, and the objective function of Equation (20) was used. Considering the efficiency and calculation accuracy, the series order, k, in the objective function was set to five. Subsequently, the temperature-independent parameters (E 0 , E ∞ , α, and β) were recalculated using the GA. The lower limit of the optimization interval was [Ê 0Ê∞αβ ](1 − ∆ 2 ), and the upper limit was [Ê 0Ê∞αβ ](1 + ∆ 2 ). The ∆ 1 and ∆ 2 were used to compensate for the rounding and fitting errors [10]. The obtained E 0 , E ∞ , α , and β are presented in Table 9. Table 9. Optimized values based on the GA. Using the initial optimized values obtained with the GA, a secondary optimization was performed using the L-M algorithm. The optimized values are shown in Table 10. Based on the optimized temperature-independent parameters, as shown in Table 10, the relaxation times at different temperatures (τ i ) were calculated, and the results are presented in Table 11. A reference temperature T r = 30°C was adopted. Based on the optimized τ i , the shift factors α T (T) at each temperature and the reduced frequency ω red could be calculated with Equations (21) and (22), respectively. Consequently, the experimental E and E data at different temperatures could be plotted at the reference temperature. Additionally, the H-N model can be employed to construct master curves based on the optimized parameters. Figure 6a,b illustrate these master curves. Furthermore, the parameters obtained using the GA and the L-M algorithm were back-calculated and plotted on the same coordinate system as the experimental data to form a Cole-Cole curve. A comparison was made to evaluate the fitting effect, as shown in Figure 6c.
τ is an internal time index that characterizes the relaxation phenomena of materials and reflects their viscoelastic nature when subjected to external forces of similar time magnitude. Temperature has a significant effect on τ, thereby affecting the viscoelastic behavior of polymer materials. The Cole-Cole curve, which represents the relationship between the storage modulus and the loss modulus, is independent of temperature and unaffected by τ. Nevertheless, its evaluation is crucial to capturing the viscoelastic characteristics of the material. Figure 6a,b illustrate the storage modulus and loss modulus over a wide frequency range. It can be observed that the storage modulus (E ) exhibited a gradual increase followed by a rapid increase before reaching a plateau as the frequency increased. Conversely, the loss modulus (E ) initially increased and then decreased as the frequency increased. During the optimization process, the relaxation time (τ i ) was set to 1 (i.e., treating each temperature as its own reference temperature), and separately fitting the parameters of E 0 , E ∞ , α, and β would result in fitting within a narrow frequency range, known as the frequency locking phenomenon. However, this narrow range fails to capture the complete viscoelastic properties of materials. To address this, it becomes necessary to shift the test data at different temperatures to construct master curves that span a broader frequency range. Hence, when fitting the relaxation time τ i at a reference temperature using the experimental data, it is necessary to calculate the shift factors α T (T) at each test temperature to the reference temperature. Figure 6c demonstrates that the accuracy of the results obtained using the L-M algorithm for the secondary fitting was slightly lower than that of the GA. This discrepancy could be attributed to the overfitting problem resulting from a large amount of data or the greater influence of initial values on the L-M algorithm. Inaccurate initial values may lead to significant deviations between the fitting results of the L-M algorithm and the actual data. Therefore, the advantage of the GA is highlighted, as it can determine the global optimal solution without relying on initial values. However, the results obtained using the GA are strongly influenced by the upper and lower bounds. Consequently, in this study, the shape characteristics and peak positions of the Cole-Cole curve were considered when setting the upper and lower bounds for the four parameters using the GA to ensure minimal error in the calculation results.
A reference temperature = 30 °C was adopted. Based on the optimized i τ , th shift factors ( ) at each temperature and the reduced frequency could be calcu lated with Equations (21) and (22), respectively. Consequently, the experimental E′ and E ′ ′ data at different temperatures could be plotted at the reference temperature. Addi tionally, the H-N model can be employed to construct master curves based on the opti mized parameters. Figure 6a,b illustrate these master curves. Furthermore, the parameter obtained using the genetic algorithm and the L-M algorithm were back-calculated and plotted on the same coordinate system as the experimental data to form a Cole-Col curve. A comparison was made to evaluate the fitting effect, as shown in Figure 6c. During the H-N model fitting process, this study initially employed the Cole-Cole curve to estimate the values of the four parameters and set the upper and lower bounds for the subsequent inversion using the GA. This approach facilitates the solution process and effectively reduces errors caused by experimental data variations, given the large amount of data involved.

Sensitivity Analysis of the Input Variables
A sensitivity analysis of the input variables of the two viscoelastic models, the 2S2P1D model and H-N model, was conducted, and a Pareto diagram was drawn, as shown in Figure 7. It can be found that the first-order sensitivity of G g and E ∞ is the largest and far exceeds that of other variables, indicating that these two input variables have a greater influence on the objective function. Therefore, it is more necessary to set the upper and lower bounds of the input variables reasonably, which is helpful for reducing unnecessary computations, improving the efficiency of the optimization algorithm, thus searching for and finding the optimal solution more effectively.

Findings and Conclusions
This paper presented the development of optimization methods for the parameters of two widely used viscoelastic models: 2P2S1D and H-N models. The genetic algorithm and L-M algorithm were utilized for parameter optimization. Based on the results of this study, the following findings can be summarized: 1. Two optimization methods were employed for the 2S2P1D model: one independent of the WLF equation utilizing the genetic algorithm and one dependent on the WLF equation utilizing the L-M algorithm. A comparison of the fitting results revealed that the optimization method independent of the WLF equation exhibited higher accuracy in parameter fitting, while the method dependent on the WLF equation demonstrated slightly lower accuracy. 2. For the 2S2P1D model, the fitting accuracy of the parameters optimized using the genetic algorithm could be further improved by secondary optimization using the L-M algorithm. This suggests that a combination of the genetic algorithm and L-M algorithm improves the overall fitting accuracy. 3. The HN model often encounters complex fractional-order derivatives during the genetic algorithm fitting, resulting in a loss of physical meaning in the solution domain and optimization failure. To address this issue, the Cole-Cole curve was initially fitted with a third-order polynomial function, followed by further fitting using the genetic algorithm with given upper and lower bounds. This improved method can increase the accuracy of the calculated results. However, the optimized values are significantly influenced by the experimental data, and it is challenging to improve the accuracy using the L-M algorithm. This may be closely related to the method used to collect the fitting data values.

Findings and Conclusions
This paper presented the development of optimization methods for the parameters of two widely used viscoelastic models: 2P2S1D and H-N models. The GA and L-M algorithm were utilized for parameter optimization. Based on the results of this study, the following findings can be summarized:

1.
Two optimization methods were employed for the 2S2P1D model: one independent of the WLF equation utilizing the GA and one dependent on the WLF equation utilizing the L-M algorithm. A comparison of the fitting results revealed that the optimization method independent of the WLF equation exhibited higher accuracy in parameter fitting, while the method dependent on the WLF equation demonstrated slightly lower accuracy.

2.
For the 2S2P1D model, the fitting accuracy of the parameters optimized using the GA could be further improved by secondary optimization using the L-M algorithm. This suggests that a combination of the GA and L-M algorithm improves the overall fitting accuracy.

3.
The H-N model often encounters complex fractional-order derivatives during the GA fitting, resulting in a loss of physical meaning in the solution domain and optimization failure. To address this issue, the Cole-Cole curve was initially fitted with a thirdorder polynomial function, followed by further fitting using the GA with given upper and lower bounds. This improved method can increase the accuracy of the calculated results. However, the optimized values are significantly influenced by the experimental data, and it is challenging to improve the accuracy using the L-M algorithm. This may be closely related to the method used to collect the fitting data values. Institutional Review Board Statement: The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Data Availability Statement: Data will be made available on request.