Different Nonlinear Regression Techniques and Sensitivity Analysis as Tools to Optimize Oil Viscosity Modeling

Four nonlinear regression techniques were explored to model gas oil viscosity on the base of Walther’s empirical equation. With the initial database of 41 primary and secondary vacuum gas oils, four models were developed with a comparable accuracy of viscosity calculation. The Akaike information criterion and Bayesian information criterion selected the least square relative errors (LSRE) model as the best one. The sensitivity analysis with respect to the given data also revealed that the LSRE model is the most stable one with the lowest values of standard deviations of derivatives. Verification of the gas oil viscosity prediction ability was carried out with another set of 43 gas oils showing remarkably better accuracy with the LSRE model. The LSRE was also found to predict better viscosity for the 43 test gas oils relative to the Aboul Seoud and Moharam model and the Kotzakoulakis and George.


Introduction
The modeling of characteristics of petroleum and its derivatives has been a subject of numerous studies [1,2]. Different regression techniques [3][4][5][6][7][8][9][10][11][12][13][14] and artificial intelligence [15,16] (machine learning, neural network) approaches have been applied to model petroleum characteristics. Nonlinear regression has been the most used approach for model parameter estimation [17]. Typically, it minimizes an objective function based on the sum of squares of errors between experimental and calculated values [17]. Usually, the models have various parameters to be determined, and sometimes multiple solutions of the objective function can be obtained. The optimal solution depends mostly on the initial guess of parameters [17]. The appropriate parameter estimation has been reported to assure by application of sensitivity analysis on the calculated parameter values [17]. The sensitivity analysis (SA) is the study of how the variation in the output of a model (numerical or otherwise) can be apportioned, qualitatively or quantitatively, to different sources of variation, and how the given model depends on the information fed into it [18]. Good modeling practice requires that the modelers provide an evaluation of the confidence in the model, possibly assessing the uncertainties associated with the modeling process and with the outcome of the model itself [18]. Originally, SA was created to deal simply with uncertainties in the input variables and model parameters. Over the course of time, the ideas have been extended to incorporate model conceptual uncertainty, that is, uncertainty in model structures, assumptions, and specifications [18]. In our recent research [14], we developed an empirical model to predict the viscosity of secondary vacuum gas oils (VGOs) that outperformed the existent empirical models published in the literature. This model was developed based on data for 24 VGOs, extending the model of Aboul-Seoud and Moharam [1] by separating the influence of the specific gravity and the average boiling point on the VGO viscosity, adopting the idea of Kotzakoulakis and George [7]. The model was validated with data for 10 additional VGOs not included in the initial database of 24 VGOs showing a better prediction ability than the model of Aboul-Seoud and Moharam [14]. In that study [14], we applied nonlinear regression using the classical approach for estimation of model parameters by minimization of the sum of squares of errors between experimental and calculated values. The viscosity measurement, however, is associated with a relatively high error (about 5% repeatability, and about 15% reproducibility) [19]. The error in viscosity measurement in our recent study [14] was found to linearly increase with the temperature of the measurement decreasing (between 5.5 and 57.8% for the temperature range 60-100 °C, being the lowest at the highest temperature).
The model parameters can be estimated not only by minimization of the sum of squares of errors between experimental and calculated values but also by minimization of the sum of absolute errors, and by minimization of the sum of relative errors [20]. Which of these nonlinear regression methods gives the best prediction is a question that needs to be investigated. For that reason, we employed data of 41 VGOs from primary and secondary origin to examine the application of four nonlinear regression methods: classical least square method, minimization of the sum of absolute errors, minimization of the sum of the squares of relative errors, and the minimization of the sum of the absolute relative errors for modeling of VGO viscosity prediction with the aim to answer the question which nonlinear regression method provides the most appropriate prediction of viscosity of VGO and other oils.
Hernández et al. [3], Hosseinifar, and Jamshidi [4], Samano et al. [17], and Alcazar, and Ancheyta [21], after the application of nonlinear regression, employed sensitivity analysis to find the most appropriate values of the model parameters. This approach was also adopted in this work and extended not only to the model parameters but also to the given data. In the works mentioned above [3,4,17,21] no sensitivity analysis with respect to the given data has been carried out.
The aim of this research is to evaluate which nonlinear regression technique is best suited to model oil viscosity and how the application of sensitivity analysis with respect to obtained model parameters and with respect to given data can assist in the selection of the most appropriate model.

Experimental Materials and Methods
Kinematic viscosity at 80 °C, specific gravity, average boiling point, refractive index, molecular weight, and aromatic ring index of the 43 VGOs from primary and secondary origin were used to develop the empirical model for prediction of viscosity applying the four nonlinear regression methods are presented in Table 1 [23].

Models
Walther's equation [24] was used as a basis for the empirical modeling of the viscosity of oils [7,10]. Mehrotra [10] proposed a correlation that has the form: Aboul-Seoud and Moharam [1] modified Equations (3)-(5) by including in it oil specific gravity and the empirical model then took the form: where 1 = 4.3414( ) 0.2 + 6.6913 and 2 = −3.7.
To estimate the components of parameter a we used four optimization methods: Method 1: Classical least squares method: Method 2: Minimization of the sum of absolute values: Method 3: Minimization the sum of squared relative errors: : Method 4: Minimization the sum of absolute relative errors:

Computational Minimization
In many cases, there are well-known specialized algorithms for global optimization. Such a case is when f is a monotone function, see for example [25,26]. On the other hand, there are many examples when the sum of squares can have several local minima, see for example [26] and references therein. In our case, we did not have any conditions guar-anteeing the convergence of an iterative process to the global extremum. In this study, one of the goals was to examine that the above-stated four methods are adequate mathematical models capable of satisfactorily describing the data. In order to do this, a minimum of the difference between measured and predicted oil viscosity (different for each model) was searched and sensitivity of model parameters with respect to the data was performed.
As an initial guess, the following modification of Aboul-Seoud and Moharam correction to Walther's model was used: If one starts the computations with the above-mentioned initial conditions, that is, when 1 = 4.3414 and 1 = −15.01620372 many overflow warnings/errors are obtained.
Using a set of quasi-random points in a five-dimensional parametric space in the neighborhood of the initial guess and calculating the values of corresponding criterion function , the computations started (for method 1) with the initial condition: More precisely, Halton sequences of quasi-random numbers, with base 2-6 were used to cover the hypercube neighborhood of the initial guess and the lengths of vertices 2. As examples, Halton squares of 20 × 20 points with bases 2, 3, and 4, 5 are plotted on Figure 1. One may compute the initial condition using initial guess and Halton points with indices and bases (0, 2), (8,3), (10,4), (3,5), and (5, 6), respectively. The discovery strategy for initial conditions of the other three methods was the same.
All computations were performed by the use of CAS Maple and NLPSolve with Modified Newton Iterative Method starting from the corresponding initial condition. The stop-criteria is the absolute difference of two consecutive iterations to be less or equal to 0.01.
Moreover, taking in mind the above fact, the appropriateness of the estimated parameters was verified by a sensitivity analysis using perturbations of model parameters in the range of ±20%, similarly as it is described by the authors of [3,4,17,21].
The gradients of Lagrangians of stated methods are calculated. Calculating the arithmetic mean and deviation of derivatives with respect to xi (for example), the standardized deviation of derivatives are interpreted as sensitivity coefficients with respect to .

Sensitivity Analysis of Least Squares Method
The classical least square problem is equivalent to the following Lagrange problem subject to The Lagrangian function for the least square method (19), (20) is Therefore, the sensitivities with respect to are be the arithmetic mean of derivatives. Let be the variance of derivatives (here the Bessel's correction is used). Standardizing, the sensitivity coefficients with respect to are obtained.
Similarly, the sensitivities with respect to and are Using both equalities in (26) and (27), is derived From (26), arithmetic mean, and variance , the sensitivity coefficients with respect to are calculated Analogously, using (28), calculated values of 1 ( ) corresponding arithmetic mean and variance 2 , the sensitivity coefficients with respect to are obtained.
Let us mark, that sometimes it is suitable to have the expressions for derivatives of L1 in terms of Lagrange multipliers : It is straightforward

Sensitivity Analysis of Absolute Value Minimization Problem
Analogously, it is suitable to consider the following constrained analog to absolute value minimization problem in method 2: The Lagrangian for the problem (33)-(36) where ji are the Lagrange multipliers, It follows from a well-known lemma from the proof of Karush-Kuhn-Tucker conditions (in fact Fritz John conditions), see [27] that if 0 is an optimal solution of the problem (33)-(36), then there exist multipliers 0 0 , 0 such that 0 0 , 0 ≥ 0, = 1, 2, 3, = 1, . . . , , not all zero, and where: where is the i-th unit vector and = (1, . . . , 1) . Let us note: Hence, one may construct a non-negative solution of the linear system setting:

Sensitivity Analysis of Squared Relative Errors
Analogously to previous subsections, the minimization problem in Method 3 is equivalent to the following Lagrange problem. The equivalent Lagrange problem is The Lagrangian is The first derivatives are (here the already calculated derivatives of 1 are used) where = 1, . . . , .
The formulas for sensitivity analysis of the sum of absolute relative errors are omitted because they are similar to the explanation in Section 2.2.7.

Results
The data in Table 1 indicate that the selected vacuum gas oils (VGO) differentiate significantly in their properties. The most important for modeling viscosity oil properties: specific gravity, and average boiling point [14] varied in the range 0.838 ÷ 1.177 for specific gravity, and 309 ÷ 488 °C for average boiling point. The VGO viscosity at 80 °C varied between 3.6 and 312.8 mm 2 /s. The Bayesian approach was used over several classical distributions to find the distribution functions of specific gravity (SG) and average boiling point (ABP). Using the Bayesian information criterion, one may conclude that the best distribution for SG is the normal distribution with mean and variance 0.98712, 0.0771899, respectively. The second and third candidates for continuous probability distribution are Gamma distribution and LogNormal distribution. The histogram and PDF (Probability Density Function) of SG data are plotted on Figure 3.  Table 1.
For the second data-ABP, using similar arguments, we again obtained the normal distribution with mean and variance 416.284, 39.0181, respectively. The histogram and PDF function are plotted in Figure 4.  Table 1. Table 2 presents data about regression coefficients for the four methods obtained after application of the Newton iterative procedure, and after the sensitivity analysis with respect to obtained parameters. These data show that the performed sensitivity analysis with respect to the model parameters in most cases led to a modification of the values of the regression coefficients. Table 3 indicates data of calculated viscosity of the 41 VGOs from Table 1, error, absolute relative error, and average abs. rel. error (AARE), also known as %AAD (average absolute deviation) by the use of the optimized values of the regression coefficients from Table 2 (model parameters after sensitivity analysis). The errors, and the %AAD were computed as shown in Equations (53) and (54) respectively: Considering the %AAD as a criterion for classification of the four studied methods the method %AAD increases in the order: Method 3 < Method 4 < Method 1 < Method 2. Table 4 shows the standardized sensitivities for the four studied estimation methods. The data in Table 4 display that Method 1 VGO under numbers 6,9,15,20,24,27 exhibited a highsensitivity for the data of ABP (x). The VGO under number 24 demonstrated high sensitivity for the data of viscosity (z). The VGOs under numbers 6, 9, 24, 27 indicate a high sensitivity for the data of SG (y). For Method 2 only the data of VGO under number 19 exhibits a high sensitivity for the data of ABP (x) and SG (y). For Method 3, the data of VGO under numbers 35 and 37 demonstrate a high sensitivity for the data of viscosity (z). The VGOs under numbers 11, 20, and 37 show high sensitivity for the data of ABP (x). The data for SG (y) of VGOs under numbers 20, and 37 display a high sensitivity. Method 4 indicates a high sensitivity for VGOs under numbers 10 and 35 for the data of viscosity (z), under number 19 for ABP (x) and SG (y). Table 5 presents data about the means and standard deviations of derivatives. It is evident from these data that Method 3 is characterized with the lowest standard deviation of derivatives followed by Method 4. Methods 1 and 2 have two and three orders of magnitude higher standard deviation of derivatives than those of Methods 3 and 4. Table 6 presents independent data (kin. viscosity at 80 °C; ABP, and SG)) for 43 gas oils to verify the capability of the four methods to predict viscosity. These data include gas oils ranging from light gas oil to VGO. The SG and ABP for this independent data set vary between 0.805 and 1.006, and between 205 and 463 °C, respectively. The kinematic viscosity varies between 0.8 and 28.1 mm2/s. The % AAD increases in the order Method 3 (18.2%) < Method 4 (28.3%) < Method 1 (61.8%) < Method 2 (67.8%).

Evaluation of the Accuracy of Viscosity Estimation by the Studied four Methods
Besides the error (53), and %AAD (54) the following additional statistical parameters were used to evaluate the accuracy of viscosity estimation by the studied four methods for the data set of Table 1 Table 7 summarizes the statistical analyses for the four studied methods employing the data in Table 1. According to the statistical parameters standard error, relative standard error Methods 1 and 2 surpass in the accuracy of viscosity prediction Methods 3 and 4. However, regarding the statistical parameters relative error, the sum of square errors, %AAD, Method 3 seems to be the best. It is difficult to distinguish the best method on the basis of the statistical parameters estimated by Equations (53)-(59). The Akaike information criterion (AIC) and Bayesian information criterion (BIC) were found capable of estimating the relative quality of a statistical method, and thus being able of providing means for model selection [13,28,29] when several models are available. Below the estimation of AIC and BIC for the four studied methods is summarized: Akaike Information Criterion. Consider the obtained errors { 1 , … , } as independent random samples from a density function ( _ | ), n = 41. Supposing normal distribution of errors: Then by the definition of likelihood function: The function L has maximum, if Again: the model with the lowest BIC score is preferred. On the base of AIC, and BIC one may conclude that Method 3 is the model with the highest quality.

Sensitivity Analysis with Respect to Given Data
Tables 4 and 5 summarizes the means and standard deviations of derivatives of the four investigated methods.
The variances 2 in the datasets of derivatives (especially with respect to yi, i = 1, ..., n, we have 2 ≈ 1408 850, respectively) are huge in the first two methods.
In Table 4 the extreme values of sensitive coefficients are marked in bold. In fact, the extreme value of the deviation of the derivative is a sign of a possible problem: the model is not suitable for certain data or the given point does not correspond to the model, etc. On the other hand, different objective functions and corresponding analyses produce different extremal values in the set of all sensitivities. Therefore, it is a good idea to perform sensitivity analysis through different objective/target functions to one and the same mathematical model and to analyze the obtained values in order to improve the model or to exclude an initially given data. In Figure 5 the distributions of sensitivity coefficients in the four methods are presented. Figure 5. The boxed-graphs of sensitivity coefficients with respect to z, x, and y for all four methods. Each boxed-graph is based on the five-number statistical characteristics; minimum, first quartile, median, third quartile, and maximum (the central rectangle represents from the first quartile to the third quartile; the segment inside the rectangle is the median; the dot is the mean). The two mentioned sets of derivatives, mentioned above, are spread out from their average value-the mean. A situation like this is possible if we did not find the extremum, or if the model function is not adequate, or the derivatives are huge in "any" small neighborhood of the extremum. In any case, the calculated values of variance are reasons to doubt the first two methods. Contrary to the third and fourth method the variances are not so huge numbers. As example we present toon Figure 6 the histogram of derivative values, computed for the first and fourth methods. Based on the arguments above, one may consider Method 3 or Method 4. Preferably, Method 3 taking in account the mean absolute percentage error.

Verification of the Viscosity Prediction Ability of the Four Studied Methods
The 43 gas oils from Table 6 were selected in such a way to cover the whole possible diversity of properties of gas oils from primary and secondary origin which can encounter in any refinery all over the world. As was already mention in section "Results" Method 3 surpassed all other methods concerning the accuracy of viscosity prediction. Table 8 summarizes the statistical analyses for the four studied methods employing the data in Table 6. These data indubitably reveal the superiority of Method 3 as the best method to model gas oil viscosity. As a supplement the oil viscosity models of Aboul Seoud and Moharam [1] (Equation (6)), and Kotzakoulakis, and George [7] (Equation (65)), which are based on Walther's equation, were verified to predict viscosity of the 43 gas oils from Table 6. They predict the 43 gas oil viscosities with % AAD of 21.8 %, and 89 % respectively proving the superiority of Method 3 model.  The model obtained by Method 3 is currently used not only to predict viscosity of gas oils but also as a tool for verification of the correctness of viscosity measurement of gas oils in LUKOIL Neftohim Burgas Research laboratory. Several times it proved its usefulness as an indicator for incorrect viscosity measurement especially when H-Oil gas oils which contain both high amount of aromatic compounds and relatively high content of waxes that makes problematic their viscosity measurement. Once HVGO viscosity at 80 °C was measured equal to 72 mm 2 /s while the model based on Method 3 reported the value of 54 mm 2 /s. The repetition of the viscosity measurement reported the value of 54 mm 2 /s.

Conclusions
The gas oil properties average boiling point and specific gravity along with modified Walther's equation and nonlinear regression techniques can be used to model oil physical property viscosity. The four nonlinear regression techniques: least squares of absolute errors, least absolute errors, least squares of relative errors, and least absolute relative errors can model gas oil viscosity. The developed gas oil viscosity models by use of the four nonlinear regression methods showed comparable accuracy of viscosity calculation for the initial base of 41 vacuum gas oils. The statistical parameters relative error, standard error, relative standard error, sum of square errors, % average absolute deviation, coefficient of determination were not in position to unequivocally select the best model. Both AIC, BIC and the standard deviations of derivatives unambiguously indicated that the model developed by nonlinear regression least squares of relative errors was the best one. The sensitivity analysis with respect to given data also revealed that the LSRE model is the most stable one with the lowest values of standard deviations of derivatives.
The LSRE model demonstrated the highest accuracy of viscosity prediction of 43 gas oils not included in the initial data base. It was also superior in oil viscosity prediction relative to other published models based on modified Walther's equation. The LSRE can be used not only to predict gas oil viscosity but also to examine the correctness of the oil viscosity measurement.

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