# Parametric Analysis of a Heavy Metal Sorption Isotherm Based on Fractional Calculus

Authors: Enrico M. Gomes; Rosana R. L. Araújo; Marcelo K. Lenzi; Fernanda R. G. B. Silva; Ervin K. Lenzi
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Heavy metals are widely recognized as being hazardous to human health and environmentally aggressive. The literature reports different approaches for lead removal, for example, water hyacinths. Heavy metal sorption isotherm modeling represents an important tool towards the study of equilibrium conditions. Fractional calculus represents a novel approach and a growing research field for process modeling, based on derivatives of arbitrary order. Recently, a novel isotherm based on fractional calculus was proposed for lead sorption using water hyacinth (Eichhornia crassipes). This paper reports a general procedure on error analysis and its influence on parameter estimation. It was applied to mathematical models based on fractional differential equations, focusing on a heavy metal novel isotherm sorption model. Parameter variance was calculated by using two different approaches (with the complete Hessian matrix and with a simplified Hessian matrix), and joint parameter confidence regions were generated, being successfully able to show that the fractional nature of the model is statistically valid.


---

## Body

## 1. Introduction

Heavy metals are widely recognized as being hazardous to human health and environmentally aggressive, being continuously generated by different chemical plants. The use of lead in the battery industry [1] is an important example. The literature reports different approaches for heavy metals removal, such as chemical precipitation [2], ion exchange [3], and electrochemical [4] and water hyacinths [5]. Mathematical models represent an essential tool for in-depth process studies, design, optimization, and control [6]. Therefore, heavy metal sorption isotherm modeling represents an important way towards the study of equilibrium conditions, which play a key role in sorption process design. The most common approach for this task consists in the use of classical models [7], such as Langmuir, Freundlich, and Redlich-Peterson among others, followed by proper parameter estimation and model discrimination analysis.Fractional calculus represents a novel approach and a growing research field for process modeling, being based on derivatives of arbitrary order [8–14]. The literature reports a broad range of applications, concerning systems engineering [15], diffusion processes [16], heat transfer [17], solid mixing [18], biological systems [19], and fluid mechanics [20] among others [21]. Recently, dos Santos et al. [22] proposed a novel isotherm based on fractional calculus for lead sorption using water hyacinths (Eichhornia crassipes). The reported isotherm can successfully predict equilibrium concentrations of lead between the aqueous solution and the water hyacinth after. The model was validated using synthetic effluent [1]. It is important to highlight that the proposed model also leads to better performances when compared to classical models (Langmuir, Freundlich, and Redlich-Peterson), which were used for sake of comparison.Error analysis represents a crucial step in model validation and further applications [23]. Recently, Joshi et al. [24] presented a detailed model analysis concerning classical sorption models. Regarding fractional-calculus-based models, Gabano and Poinot [25], Khemane et al. [26], and Isfer et al. [15] report the calculation of parametric variance.It is important to state that one may identify a mathematical model of fractional order for a given set of experimental data [27]; however, only a precise error analysis can ensure that the derivative is in fact fractional. If the variance of the estimated fractional order of the derivative is large enough, the fractional order can be statistically regarded as integer order for a given confidence level. Consequently, the analysis of parameter joint confidence region becomes an essential tool, as the region indicates, for a given confidence level, the possible set of parameters that could generate the experimental data [28]. To the best of our knowledge, the generation and analysis of joint confidence regions have not yet been reported for models based on fractional calculus. This paper reports a detailed study on an error analysis procedure applied to mathematical models based on fractional differential equations. After the development of a theoretical framework concerning parameter estimation, parameter variance estimation and joint confidence region determination, the fractional model proposed by dos Santos [22] was used as a case study for validation purposes.


## 2. Theoretical Framework

### 2.1. Mathematical Model

Further details regarding the experimental data set can be obtained from dos Santos and Lenzi [1]. It needs to be highlighted that experimental data was normalized in the interval [0,1] for proper parameter estimation [22, 29]. The mathematical model used in this work was firstly proposed by dos Santos et al. [22] for describing the lead equilibrium sorption. According to the authors, a large number of experimental results on equilibrium systems dealing with heavy metals have the following behavior: when the heavy metal concentration in the aqueous phase is low, the equilibrium concentration in the solid phase may largely change for a given modification in the concentration of the fluid phase. On the other hand, for higher concentrations of the heavy metal in the fluid phase, the equilibrium concentration in the solid may be less sensitive, indicating some kind of saturation. These features resemble to a certain degree in an exponential behavior, obtained, for example, from first-order differential equations. Consequently, an exponential model for heavy metal sorption isotherm, as given by (1), can explain some normalized experimental results, where parameters θ1 and θ2 depend on the sorption process features, like the type of heavy metal, solid matrix used for sorption process, and so on:
(1)θ1·dydx+θ2·y=1,y(x=0)=0⟹y=(1θ2)·(1-e-(θ2/θ1)·x).Therefore, by using fractional calculus, the previous model can be generalized to (2), by considering a fractional order θ3 for the differential equation. One can note that parameter θ3 also depends on the features of the sorption process:(2)θ1·dθ3qdxθ3+θ2·y=1,y(x=0)=0⟹y=(1θ1)·ε0(x,-θ2θ1;θ3,θ3+1),Epsilonfunction:ε0(m1,m2;m3,m4)=(m1)(m4-1)·Em3,m4(m2·(m1)m3),Mittag-Lefflerfunction:En1,n2(n3)=∑j=0∞(n3)jΓ(n1·j+n2),
where ε0 is the 0th-order Epsilon function defined by Podlubny [30], which uses the Mittag-Leffler function; mi  i={1,…,4} and ni  i={1,…,3} are dummy variables. Details regarding the Gamma function (Γ) can be found in the appendix.


### 2.2. Parameter Estimation

Parameter estimation was carried out using a genetic algorithm procedure as reported by Isfer et al. [15]. More specifically, the initial population consisted of 250 sets of values for parameters θ_={θ1;θ2;θ3}, which iterated until the difference of each parameter θi, the best set of two consecutive iterations, was lower than 10-6. Crossover and mutation probabilities were of 80% assuring a good macroscopic search and of 10% assuring a good microscopic (refinement), respectively. In order to avoid a local optimum solution, estimation was performed using different initial populations. Parameters θ1, θ2, and θ3 were estimated using (3) as the objective function of the optimization problem, representing a normalized quadratic error analysis [31]. Experimental variances were considered constant and equal to δyE2 for all experiments and experimental covariances were assumed equal to zero. Consequently, matrix [VyE__](NE×NE) is a diagonal matrix. According to Bard [32], δyE2 can be approximated by (4):
(3)FOBJ=[(yE_-yM_)T](1×NE)·[(VyE__)-1](NE×NE)·(yE_-yM_)(NE×1)=(1δyE2)·∑i=1NE(yiE-yiM(θ_;xi))2,(4)δyE2=∑i=1NE(yiE-yiM(θ_;xi))2NE-NP.


### 2.3. Parameter Variance

According to Bard [32], for the objective function defined by (3), the parametric variance matrix [Vθ__](NP×NP) is given by (5), which uses the Hessian matrix [Hθ__](NP×NP) and matrix [Gy__](NP×NE), given by (6):
(5)[Vθ__](NP×NP)=[(Hθ__)-1](NP×NP)·([Gy__](NP×NE)·[VyE__](NE×NE)·[GyT__](NE×NP))·[(Hθ__)-1](NP×NP)(6)[Hθ__](NP×NP)=[∂2FOBJ∂θ12⋯∂2FOBJ∂θ1∂θNP⋮⋱⋮∂2FOBJ∂θ1∂θNP⋯∂2FOBJ∂θNP2],[Gy__](NP×NE)=[∂2FOBJ∂θ1∂y1E⋯∂2FOBJ∂θ1∂yNEE⋮⋱⋮∂2FOBJ∂θNP∂y1E⋯∂2FOBJ∂θNP∂yNEE].It can be observed that matrix[Gy__](NP×NE) can also be obtained by using a sensitivity matrix [23], [B__](NE×NP), which is given by (8):
(7)[Gy__](NP×NE)=2·[(B__)T](NP×NE)·[(VyE__)-1](NE×NE),[(B__)T](NP×NE)=(12)·[Gy__](NP×NE)·[VyE__](NE×NE),(8)[B__](NE×NP)=[∂y1M∂θ1⋯∂y1M∂θNP⋮⋱⋯∂yNEM∂θ1⋯∂yNEM∂θNP].Also, according to Bard [32], the elements [hij] of the Hessian matrix [Hθ__](NP×NP) are given by (9)∂2FOBJ∂θi∂θj=2·[(∂yM_∂θi)T](1×NE)·[(VyE__)-1](NE×NE)·(∂yM_∂θj)(NE×1)︸TERM1-2·(∂2yM_∂θi∂θj)(1×NE)·[(VyE__)-1](NE×NE)·(yE_-yM_)(NE×1)︸TERM2.For linear estimation problems, TERM2 automatically vanishes. However, this is not the case for nonlinear problems, but, according to Alberton et al. [23], this term can be neglected in some scenarios. Therefore, the parametric variance matrix is usually approximated by (10)[Vθ__](NP×NP)=([(B__)T](NP×NE)·[(VyE__)-1](NE×NE)·[B__](NE×NP))-1=4·(([Gy__](NP×NE)·[VyE__](NE×NE)·[GyT__](NE×NP)))-1.


### 2.4. Parameter Confidence Interval and Joint Confidence Region

According to Himmelblau [33], parameter confidence intervals can be obtained by using (11) for a given confidence level of 100·(1-α)%, but the use different values of ρ has been reported. For a small number of experimental data, the use of ρ=t(1-α/2),(NE-NP), obtained from Student’s t-distribution, is recommended. On the other hand, ρ=z(1-α/2), obtained from the normal distribution, can be used for a large number of experimental data:
(11)θi#-ρ·δθi<θi#<θi#+ρ·δθi.Based on parametric variance, joint confidence region also needs to be obtained. For a given confidence level, this region provides the set of parameters that could actually generate the experimental data set. This is an important analysis tool, as although a given set of parameter may be within the confidence interval, it may not be inside the joint confidence region [33]. Usually, these joint confidence regions can be obtained by (12), which considers a linearization of the estimation problem [28]. Thus, a key issue to be addressed is concerns on the influence of the approach used to calculate matrix [Vθ__](NP×NP) ((5) or (10)) on the shape of the confidence region:
(12)[(θ_-θ_#)T](1×NP)·[(Vθ__)-1](NP×NP)·[(θ_-θ_#)](NP×1)≤FOBJ(θ_#)·NPNE-NP·FNP,(NE-NP)(1-α).It needs to be stressed that (12) is only an approximation for nonlinear problems. A much more realistic approach [28, 34] given by (13) which considers the intrinsic nonlinear features of the estimation problem, leading to joint confidence regions usually larger than the ones obtained by (12):
(13)FOBJ(θ_)≤FOBJ(θ_#)·(1+NPNE-NP·FNP,(NE-NP)1-α).


### 2.5. Model Prediction Variance

Pinto and Schwaab [35] report that the variance of the model predictions of the experimental data used for parameter estimation is given by (14). The main feature of this equation is that not only experimental and modeling error themselves are taken into account, but also the correlation between them is also considered for the variance prediction:(14)[VyM#__](NE×NE)=[B__](NE×NP)·[Vθ__](NP×NP)·[(B__)T](NP×NE)︸modelingerrors+[VyE__](NE×NE)︸experimentalerrors-([B__](NE×NP)·([(Hθ__)-1](NP×NP)·[Gy__](NP×NE)·[VyE__](NE×NE))+([(Hθ__)-1](NP×NP)·[Gy__](NP×NE)·[VyE__](NE×NE))T·[(B__)T](NP×NE))︸correlationbetweenexperimentalandmodelingerrors.For prediction of the variance either of future experiments (for a given value ofxi) or available experimental data not used during parameter estimation, (15) needs to be used. Only used data for parameter estimation contribute to the correlation between experimental and modeling errors:
(15)δyiM⊗2=([B⊗__](1×NP)·[Vθ__](NP×NP)·[(B⊗__)T](NP×1))+(δyE2),
where [B⊗__](1×NP) is the model gradient vector given by (16) and evaluated for xi:
(16)[B⊗__](1×NP)=[∂yM∂θ1∂yM∂θ2⋯∂yM∂θNP].


### 2.6. Analytical Expressions

By using the definitions ofε0 and the Mittag-Leffler function, the mathematical model given by (2) can be rewritten as
(17)yiM=∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)].The objective function used for estimating the parameter setθ_={θ1;θ2;θ3} is given by (18)
(18)FOBJ=∑i=1NE(yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))2.The elements of matrix[B__](NE×NP) and matrix [B⊗__](1×NP) are given by
(19)∂yiM∂θ1=∑j=0∞[-(-1)j·(θ2)j·xi((θ3)·(j+1))·(j+1)(θ1)j+2·Γ((θ3)·(j+1)+1)],∂yiM∂θ2=∑j=0∞[(-1)j·(θ2)(j-1)·xi((θ3)·(j+1))·j(θ1)j+1·Γ((θ3)·(j+1)+1)],∂yiM∂θ3=∑j=0∞[((θ1)j+1·Γ((θ3)·(j+1)+1))-1((-1)j·(θ2)j·xi((θ3)·(j+1))·(j+1)·(ln(xi)-Ψ((θ3)·(j+1)+1))(-1)j)×((θ1)j+1·Γ((θ3)·(j+1)+1))-1].The objective function gradient is given by(20)∂FOBJ∂θ1=∑i=1NE(-2·(yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)])))·∂yiM∂θ1,∂FOBJ∂θ2=∑i=1NE(-2·(yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)])))·∂yiM∂θ2,∂FOBJ∂θ3=∑i=1NE(-2·(yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)])))·∂yiM∂θ3.The elements of matrix[Gy__](NP×NE) are obtained by
(21)∂2FOBJ∂θ1∂yiE=-2·(∑j=0∞[-(-1)j·(θ2)j·xi((θ3)·(j+1))·(j+1)(θ1)j+2·Γ((θ3)·(j+1)+1)]),∂2FOBJ∂θ2∂yiE=-2·(∑j=0∞[(-1)j·(θ2)(j-1)·xi((θ3)·(j+1))·j(θ1)j+1·Γ((θ3)·(j+1)+1)]),∂2FOBJ∂θ3∂yiE=-2·(∑j=0∞[((θ1)j+1·Γ((θ3)·(j+1)+1))-1((-1)j·(θ2)j·xi((θ3)·(j+1))·(j+1)·(ln(xi)-Ψ((θ3)·(j+1)+1))xi((θ3)·(j+1)))×((θ1)j+1·Γ((θ3)·(j+1)+1))-1]∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))·(j+1)).Finally, the elements of[Hθ__](NP×NP) result from(22)∂2FOBJ∂θ12=∑i=1NE((2·(∂yiM∂θ1)2)-2·((yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))·∑j=0∞[×((θ1)j+3·Γ((θ3)·(j+1)+1))-1((-1)j·(θ2)j·xi((θ3)·(j+1))·((j+1)2+(j+1))xi((θ3)·(j+1)))×((θ1)j+3·Γ((θ3)·(j+1)+1))-1](∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))),∂2FOBJ∂θ22=∑i=1NE((2·(∂yiM∂θ2)2)-2·((yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))·∑j=0∞[(-1)j·(θ2)(j-2)·xi((θ3)·(j+1))·(j2-j)(θ1)j+1·Γ((θ3)·(j+1)+1)])),∂2FOBJ∂θ32=∑i=1NE((2·(∂yiM∂θ3)2)-2·((yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))·∑j=0∞[((-1)j·(θ2)j·xi((θ3)·(j+1))·(j+1)2·((ln(xi)-Ψ((θ3)·(j+1)+1))2-Ψ(1,(θ3)·(j+1)+1)((ln(xi)-Ψ((θ3)·(j+1)+1))2)xi((θ3)·(j+1)))×((θ1)j+1·Γ((θ3)·(j+1)+1))-1](∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))),∂2FOBJ∂θ1∂θ2=∑i=1NE((2·(∂yiM∂θ1)·(∂yiM∂θ2))-2·((yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))·∑j=0∞[((θ1)j+2·Γ((θ3)·(j+1)+1))-1-((-1)j·(θ2)(j-1)·xi((θ3)·(j+1))·(j·(j+1))(-1)j·(θ2)(j-1)·xi((θ3)·(j+1)))×((θ1)j+2·Γ((θ3)·(j+1)+1))-1](∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))),∂2FOBJ∂θ1∂θ3=∑i=1NE((2·(∂yiM∂θ1)·(∂yiM∂θ3))-2·((yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))·∑j=0∞[×((θ1)j+2·Γ((θ3)·(j+1)+1))-1((-1)j·(θ2)j·xi((θ3)·(j+1))·(j+1)2·(-ln(xi)+Ψ((θ3)·(j+1)+1))xi((θ3)·(j+1)))×((θ1)j+2·Γ((θ3)·(j+1)+1))-1](∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))),∂2FOBJ∂θ2∂θ3=∑i=1NE((2·(∂yiM∂θ2)·(∂yiM∂θ3))-2·((yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))·∑j=0∞[×((θ1)j+1·Γ((θ3)·(j+1)+1))-1(j)-1((-1)j·(θ2)(j-1)·xi((θ3)·(j+1))·(j·(j+1))·(ln(xi)-Ψ((θ3)·(j+1)+1))xi(n)(-1)j·(θ2)(j-1)·xi((θ3)·(j+1)))×((θ1)j+1·Γ((θ3)·(j+1)+1))-1](∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))).Further details regarding Gamma (Γ) and Psi (Ψ) functions can be found in the appendix.


## 2.1. Mathematical Model

Further details regarding the experimental data set can be obtained from dos Santos and Lenzi [1]. It needs to be highlighted that experimental data was normalized in the interval [0,1] for proper parameter estimation [22, 29]. The mathematical model used in this work was firstly proposed by dos Santos et al. [22] for describing the lead equilibrium sorption. According to the authors, a large number of experimental results on equilibrium systems dealing with heavy metals have the following behavior: when the heavy metal concentration in the aqueous phase is low, the equilibrium concentration in the solid phase may largely change for a given modification in the concentration of the fluid phase. On the other hand, for higher concentrations of the heavy metal in the fluid phase, the equilibrium concentration in the solid may be less sensitive, indicating some kind of saturation. These features resemble to a certain degree in an exponential behavior, obtained, for example, from first-order differential equations. Consequently, an exponential model for heavy metal sorption isotherm, as given by (1), can explain some normalized experimental results, where parameters θ1 and θ2 depend on the sorption process features, like the type of heavy metal, solid matrix used for sorption process, and so on:
(1)θ1·dydx+θ2·y=1,y(x=0)=0⟹y=(1θ2)·(1-e-(θ2/θ1)·x).Therefore, by using fractional calculus, the previous model can be generalized to (2), by considering a fractional order θ3 for the differential equation. One can note that parameter θ3 also depends on the features of the sorption process:(2)θ1·dθ3qdxθ3+θ2·y=1,y(x=0)=0⟹y=(1θ1)·ε0(x,-θ2θ1;θ3,θ3+1),Epsilonfunction:ε0(m1,m2;m3,m4)=(m1)(m4-1)·Em3,m4(m2·(m1)m3),Mittag-Lefflerfunction:En1,n2(n3)=∑j=0∞(n3)jΓ(n1·j+n2),
where ε0 is the 0th-order Epsilon function defined by Podlubny [30], which uses the Mittag-Leffler function; mi  i={1,…,4} and ni  i={1,…,3} are dummy variables. Details regarding the Gamma function (Γ) can be found in the appendix.


## 2.2. Parameter Estimation

Parameter estimation was carried out using a genetic algorithm procedure as reported by Isfer et al. [15]. More specifically, the initial population consisted of 250 sets of values for parameters θ_={θ1;θ2;θ3}, which iterated until the difference of each parameter θi, the best set of two consecutive iterations, was lower than 10-6. Crossover and mutation probabilities were of 80% assuring a good macroscopic search and of 10% assuring a good microscopic (refinement), respectively. In order to avoid a local optimum solution, estimation was performed using different initial populations. Parameters θ1, θ2, and θ3 were estimated using (3) as the objective function of the optimization problem, representing a normalized quadratic error analysis [31]. Experimental variances were considered constant and equal to δyE2 for all experiments and experimental covariances were assumed equal to zero. Consequently, matrix [VyE__](NE×NE) is a diagonal matrix. According to Bard [32], δyE2 can be approximated by (4):
(3)FOBJ=[(yE_-yM_)T](1×NE)·[(VyE__)-1](NE×NE)·(yE_-yM_)(NE×1)=(1δyE2)·∑i=1NE(yiE-yiM(θ_;xi))2,(4)δyE2=∑i=1NE(yiE-yiM(θ_;xi))2NE-NP.


## 2.3. Parameter Variance

According to Bard [32], for the objective function defined by (3), the parametric variance matrix [Vθ__](NP×NP) is given by (5), which uses the Hessian matrix [Hθ__](NP×NP) and matrix [Gy__](NP×NE), given by (6):
(5)[Vθ__](NP×NP)=[(Hθ__)-1](NP×NP)·([Gy__](NP×NE)·[VyE__](NE×NE)·[GyT__](NE×NP))·[(Hθ__)-1](NP×NP)(6)[Hθ__](NP×NP)=[∂2FOBJ∂θ12⋯∂2FOBJ∂θ1∂θNP⋮⋱⋮∂2FOBJ∂θ1∂θNP⋯∂2FOBJ∂θNP2],[Gy__](NP×NE)=[∂2FOBJ∂θ1∂y1E⋯∂2FOBJ∂θ1∂yNEE⋮⋱⋮∂2FOBJ∂θNP∂y1E⋯∂2FOBJ∂θNP∂yNEE].It can be observed that matrix[Gy__](NP×NE) can also be obtained by using a sensitivity matrix [23], [B__](NE×NP), which is given by (8):
(7)[Gy__](NP×NE)=2·[(B__)T](NP×NE)·[(VyE__)-1](NE×NE),[(B__)T](NP×NE)=(12)·[Gy__](NP×NE)·[VyE__](NE×NE),(8)[B__](NE×NP)=[∂y1M∂θ1⋯∂y1M∂θNP⋮⋱⋯∂yNEM∂θ1⋯∂yNEM∂θNP].Also, according to Bard [32], the elements [hij] of the Hessian matrix [Hθ__](NP×NP) are given by (9)∂2FOBJ∂θi∂θj=2·[(∂yM_∂θi)T](1×NE)·[(VyE__)-1](NE×NE)·(∂yM_∂θj)(NE×1)︸TERM1-2·(∂2yM_∂θi∂θj)(1×NE)·[(VyE__)-1](NE×NE)·(yE_-yM_)(NE×1)︸TERM2.For linear estimation problems, TERM2 automatically vanishes. However, this is not the case for nonlinear problems, but, according to Alberton et al. [23], this term can be neglected in some scenarios. Therefore, the parametric variance matrix is usually approximated by (10)[Vθ__](NP×NP)=([(B__)T](NP×NE)·[(VyE__)-1](NE×NE)·[B__](NE×NP))-1=4·(([Gy__](NP×NE)·[VyE__](NE×NE)·[GyT__](NE×NP)))-1.


## 2.4. Parameter Confidence Interval and Joint Confidence Region

According to Himmelblau [33], parameter confidence intervals can be obtained by using (11) for a given confidence level of 100·(1-α)%, but the use different values of ρ has been reported. For a small number of experimental data, the use of ρ=t(1-α/2),(NE-NP), obtained from Student’s t-distribution, is recommended. On the other hand, ρ=z(1-α/2), obtained from the normal distribution, can be used for a large number of experimental data:
(11)θi#-ρ·δθi<θi#<θi#+ρ·δθi.Based on parametric variance, joint confidence region also needs to be obtained. For a given confidence level, this region provides the set of parameters that could actually generate the experimental data set. This is an important analysis tool, as although a given set of parameter may be within the confidence interval, it may not be inside the joint confidence region [33]. Usually, these joint confidence regions can be obtained by (12), which considers a linearization of the estimation problem [28]. Thus, a key issue to be addressed is concerns on the influence of the approach used to calculate matrix [Vθ__](NP×NP) ((5) or (10)) on the shape of the confidence region:
(12)[(θ_-θ_#)T](1×NP)·[(Vθ__)-1](NP×NP)·[(θ_-θ_#)](NP×1)≤FOBJ(θ_#)·NPNE-NP·FNP,(NE-NP)(1-α).It needs to be stressed that (12) is only an approximation for nonlinear problems. A much more realistic approach [28, 34] given by (13) which considers the intrinsic nonlinear features of the estimation problem, leading to joint confidence regions usually larger than the ones obtained by (12):
(13)FOBJ(θ_)≤FOBJ(θ_#)·(1+NPNE-NP·FNP,(NE-NP)1-α).


## 2.5. Model Prediction Variance

Pinto and Schwaab [35] report that the variance of the model predictions of the experimental data used for parameter estimation is given by (14). The main feature of this equation is that not only experimental and modeling error themselves are taken into account, but also the correlation between them is also considered for the variance prediction:(14)[VyM#__](NE×NE)=[B__](NE×NP)·[Vθ__](NP×NP)·[(B__)T](NP×NE)︸modelingerrors+[VyE__](NE×NE)︸experimentalerrors-([B__](NE×NP)·([(Hθ__)-1](NP×NP)·[Gy__](NP×NE)·[VyE__](NE×NE))+([(Hθ__)-1](NP×NP)·[Gy__](NP×NE)·[VyE__](NE×NE))T·[(B__)T](NP×NE))︸correlationbetweenexperimentalandmodelingerrors.For prediction of the variance either of future experiments (for a given value ofxi) or available experimental data not used during parameter estimation, (15) needs to be used. Only used data for parameter estimation contribute to the correlation between experimental and modeling errors:
(15)δyiM⊗2=([B⊗__](1×NP)·[Vθ__](NP×NP)·[(B⊗__)T](NP×1))+(δyE2),
where [B⊗__](1×NP) is the model gradient vector given by (16) and evaluated for xi:
(16)[B⊗__](1×NP)=[∂yM∂θ1∂yM∂θ2⋯∂yM∂θNP].


## 2.6. Analytical Expressions

By using the definitions ofε0 and the Mittag-Leffler function, the mathematical model given by (2) can be rewritten as
(17)yiM=∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)].The objective function used for estimating the parameter setθ_={θ1;θ2;θ3} is given by (18)
(18)FOBJ=∑i=1NE(yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))2.The elements of matrix[B__](NE×NP) and matrix [B⊗__](1×NP) are given by
(19)∂yiM∂θ1=∑j=0∞[-(-1)j·(θ2)j·xi((θ3)·(j+1))·(j+1)(θ1)j+2·Γ((θ3)·(j+1)+1)],∂yiM∂θ2=∑j=0∞[(-1)j·(θ2)(j-1)·xi((θ3)·(j+1))·j(θ1)j+1·Γ((θ3)·(j+1)+1)],∂yiM∂θ3=∑j=0∞[((θ1)j+1·Γ((θ3)·(j+1)+1))-1((-1)j·(θ2)j·xi((θ3)·(j+1))·(j+1)·(ln(xi)-Ψ((θ3)·(j+1)+1))(-1)j)×((θ1)j+1·Γ((θ3)·(j+1)+1))-1].The objective function gradient is given by(20)∂FOBJ∂θ1=∑i=1NE(-2·(yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)])))·∂yiM∂θ1,∂FOBJ∂θ2=∑i=1NE(-2·(yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)])))·∂yiM∂θ2,∂FOBJ∂θ3=∑i=1NE(-2·(yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)])))·∂yiM∂θ3.The elements of matrix[Gy__](NP×NE) are obtained by
(21)∂2FOBJ∂θ1∂yiE=-2·(∑j=0∞[-(-1)j·(θ2)j·xi((θ3)·(j+1))·(j+1)(θ1)j+2·Γ((θ3)·(j+1)+1)]),∂2FOBJ∂θ2∂yiE=-2·(∑j=0∞[(-1)j·(θ2)(j-1)·xi((θ3)·(j+1))·j(θ1)j+1·Γ((θ3)·(j+1)+1)]),∂2FOBJ∂θ3∂yiE=-2·(∑j=0∞[((θ1)j+1·Γ((θ3)·(j+1)+1))-1((-1)j·(θ2)j·xi((θ3)·(j+1))·(j+1)·(ln(xi)-Ψ((θ3)·(j+1)+1))xi((θ3)·(j+1)))×((θ1)j+1·Γ((θ3)·(j+1)+1))-1]∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))·(j+1)).Finally, the elements of[Hθ__](NP×NP) result from(22)∂2FOBJ∂θ12=∑i=1NE((2·(∂yiM∂θ1)2)-2·((yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))·∑j=0∞[×((θ1)j+3·Γ((θ3)·(j+1)+1))-1((-1)j·(θ2)j·xi((θ3)·(j+1))·((j+1)2+(j+1))xi((θ3)·(j+1)))×((θ1)j+3·Γ((θ3)·(j+1)+1))-1](∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))),∂2FOBJ∂θ22=∑i=1NE((2·(∂yiM∂θ2)2)-2·((yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))·∑j=0∞[(-1)j·(θ2)(j-2)·xi((θ3)·(j+1))·(j2-j)(θ1)j+1·Γ((θ3)·(j+1)+1)])),∂2FOBJ∂θ32=∑i=1NE((2·(∂yiM∂θ3)2)-2·((yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))·∑j=0∞[((-1)j·(θ2)j·xi((θ3)·(j+1))·(j+1)2·((ln(xi)-Ψ((θ3)·(j+1)+1))2-Ψ(1,(θ3)·(j+1)+1)((ln(xi)-Ψ((θ3)·(j+1)+1))2)xi((θ3)·(j+1)))×((θ1)j+1·Γ((θ3)·(j+1)+1))-1](∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))),∂2FOBJ∂θ1∂θ2=∑i=1NE((2·(∂yiM∂θ1)·(∂yiM∂θ2))-2·((yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))·∑j=0∞[((θ1)j+2·Γ((θ3)·(j+1)+1))-1-((-1)j·(θ2)(j-1)·xi((θ3)·(j+1))·(j·(j+1))(-1)j·(θ2)(j-1)·xi((θ3)·(j+1)))×((θ1)j+2·Γ((θ3)·(j+1)+1))-1](∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))),∂2FOBJ∂θ1∂θ3=∑i=1NE((2·(∂yiM∂θ1)·(∂yiM∂θ3))-2·((yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))·∑j=0∞[×((θ1)j+2·Γ((θ3)·(j+1)+1))-1((-1)j·(θ2)j·xi((θ3)·(j+1))·(j+1)2·(-ln(xi)+Ψ((θ3)·(j+1)+1))xi((θ3)·(j+1)))×((θ1)j+2·Γ((θ3)·(j+1)+1))-1](∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))),∂2FOBJ∂θ2∂θ3=∑i=1NE((2·(∂yiM∂θ2)·(∂yiM∂θ3))-2·((yiE-(∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))·∑j=0∞[×((θ1)j+1·Γ((θ3)·(j+1)+1))-1(j)-1((-1)j·(θ2)(j-1)·xi((θ3)·(j+1))·(j·(j+1))·(ln(xi)-Ψ((θ3)·(j+1)+1))xi(n)(-1)j·(θ2)(j-1)·xi((θ3)·(j+1)))×((θ1)j+1·Γ((θ3)·(j+1)+1))-1](∑j=0∞[(-1)j·(θ2)j·xi((θ3)·(j+1))(θ1)j+1·Γ((θ3)·(j+1)+1)]))).Further details regarding Gamma (Γ) and Psi (Ψ) functions can be found in the appendix.


## 3. Results

Heavy metal sorption processes involve a solid phase, which in many scenarios consist in irregular and disordered structures; consequently, many classical models usually cannot adequately explain observed experimental behavior. Towards this, fractional calculus plays a key role in mathematical modeling of transport phenomena in irregular structure, mainly due to memory effects as already reported in the literature [36, 37]. It is important to note that the model proposed by dos Santos et al. [22] presents a novel characteristic as according to the value of parameter θ3, the Epsilon function turns into a different mathematical function; for example, if θ3 is equal to 1, an exponential form is obtained [30].Table1 presents a summary of the parameter estimation procedure. The row EPSILON 1 presents the results considering (5) to calculate the parameter variance matrix. The results in row EPSILON 2 considered (10) to calculate the parameter variance matrix. Firstly, it can be seen that the model reported by dos Santos et al. [22] presents a very good correlation coefficient and a low value of the objective function.Table 1
Parameter estimation results@.
Parameter
δ
θ
i
θ
i
±
t
(
1
-
(
α
/
2
)
)
,
(
NE
-
NP
)
·
δ
θ
i  θi±z(1-(α/2))·δθi
[
V
θ
_
_
]
[
V
r
θ
_
_
]
r
F
OBJ
θ
1
0.169
2.66
·
10
-
3
0.169 ± 0.0110.169 ± 0.005
[
7.08
·
1
0
-
6
-
2.02
·
1
0
-
6
-
5.36
·
1
0
-
6
-
2.02
·
1
0
-
6
3.83
·
1
0
-
4
9.37
·
1
0
-
5
-
5.36
·
1
0
-
6
9.37
·
1
0
-
5
2.89
·
1
0
-
5
]
[
1
-
0.039
-
0.38
-
0.039
1
0.89
-
0.38
0.89
1
]
0.999
1.38
·
10
-
3
EPSILON 1
θ
2
0.884
1.96
·
10
-
2
0.884 ± 0.0840.884 ± 0.038
θ
3
0.488
5.38
·
10
-
3
0.488 ± 0.0230.488 ± 0.010
θ
1
0.169
3.12
·
10
-
2
0.169 ± 0.1340.169 ± 0.061
[
9.74
·
1
0
-
4
-
1.57
·
1
0
-
3
-
1.34
·
1
0
-
3
-
1.57
·
1
0
-
3
3.31
·
1
0
-
3
1.93
·
1
0
-
3
-
1.34
·
1
0
-
3
1.93
·
1
0
-
3
2.05
·
1
0
-
3
]
[
1
-
0.87
-
0.95
-
0.87
1
0.74
-
0.95
0.74
1
]
EPSILON 2
θ
2
0.884
5.76
·
10
-
2
0.884 ± 0.2470.884 ± 0.112
θ
3
0.488
4.52
·
10
-
2
0.488 ± 0.1940.488 ± 0.088
@100 · (1 − α)% = 95%: t(1-(α/2)),(NE-NP)=t0.975;(5-3)=t0.975;2=4.3030; t(1-(α/2)),(NE-NP)=t0.975;(5-2)=t0.975;3=3.182; and t(1-(α/2)),(∞)=z(1-(α/2))=z0.975=1.96.It can also be observed that the parameter variance is small when compared to the value of the parameter; consequently, all the parameters can be considered significant. This fact occurred independently of the approach used for the calculations of the parametric variance. This conclusion is also obtained after analyzing the parameter confidence intervals calculated for 95% of confidence level, regardless of the value used forρ in (11). It can be seen that all parameters are significant as they are all statistically different from zero. However, it is important to emphasize that for this nonlinear parameter estimation study, the choice of the approach used to calculate the parameter variance led to differences of one order of magnitude for parameters θ1 and θ3. Besides, some interesting effects regarding parameter correlation showed up; more specifically, by choosing the simplified approach, correlation between parameters θ1 and θ3 and parameters θ1 and θ2 considerably increased.Figure1 presents the experimental values plotted against the model predictions. Bar errors were obtained using (4) to obtain an approximation of the experimental error and (14) to obtain the error of the model predictions of the experimental data used in the parameter estimation task. It can be observed that model predictions are in good agreement with the experimental values as the points are close to the straight line, which indicates the ideal case that model predictions are equal to the experimental values. Figure 2 presents a comparison of experimental data and model predictions plotted against the independent variable. As mentioned before, experimental error was calculated using (4). Again, one can see that the fractional model adequately describes the experimental data, especially for low concentrations of lead in the aqueous solutions, where low variations cause large variations in the lead concentration in the hyacinth. Besides, Figure 2 also presents the confidence interval limits (1·σ) considering (15) to calculate the model prediction error and using the complete Hessian to calculate [Vθ__](NP×NP).Figure 1
Experimental values versus model predictions.Figure 2
Isotherm behavior.Parameterθ3 plays a key role in the model because, as mentioned before, according to its value the equation can assume a different form. Besides, it is the order of the fractional differential equation (see (2)). Therefore, it is important to prove that θ3 is statistically different from 1, otherwise an integer order differential equation would be a model with the same efficiency. Initially, this analysis considers the confidence interval presented in Table 1. However, this may not be enough, as the parameter joint confidence region plays a key role. Figure 3 presents the region obtained by (12) considering the parametric variance matrix obtained by (5) (solid region) and considering the parametric variance matrix obtained by (9) (wired region). One can observe that θ3 is different from 1; consequently, the use of fractional order derivative is significant. Secondly, the approach used for calculating the parameter variance significantly influences the parameter confidence region as it can be seen by the difference in size and shape of the ellipsoids. This is an important consequence as a set of parameters which is inside the wired region may be outside the solid region; therefore, it may not be statistically significant.Figure 3
Confidence region—Approach 1.Finally, it needs to be remembered that the joint confidence regions shown in Figure3 were obtained by a simplified though useful approach. For nonlinear parameter estimation problems, a more accurate parameter confidence region is obtained by (13), which is presented by Figure 4. It is essential to stress that the size and the shape of the region may considerably change as (13) preserves the nonlinear features of the problem, remembering that (12) is somehow a linearization of (13). Moreover, it must be emphasized that parameter θ3 is still lower than 1 (see Figure 4); consequently, the model nature and experimental data behavior can be adequately described by a fractional order model.Figure 4
Confidence region—Approach 2.It is important to analyze the variance of the model predictions of future experiments, which can be obtained by (15). Figure 5 presents the variance predictions plotted against the independent variable, that is, lead concentration in the aqueous solution. The variances were calculated using the different parametric variance matrixes; that is, Epsilon1 used (5) and Epsilon2 used (10). One can observe that using (10) the variance of future experiments predictions considerably changes. It must be remembered that, although often used, (10) is an approximation of the calculation of the parametric variance. Therefore, for nonlinear problems as the one present in fractional calculus applications, (5) should be chosen. It is also important to note that the minimum values of variance are obtained for lower values of the independent variable. For higher values, the variance tends to increase.Figure 5
Variance behavior for different values of the independent variable.


## 4. Conclusions

The availability of mathematical models plays a key role in understanding heavy metal equilibrium phenomena. The literature reports different approaches for modeling the sorption process, particularly lead. Recently, a model based on fractional calculus was proposed to describe experimental data concerning lead sorption by hyacinths. This paper reported the use of an error analysis procedure for mathematical models based on fractional order differential equations in order to show that the fractional order is in fact fractional.Parametric variance matrix was calculated by two different approaches, one considering the complete Hessian matrix and the other considering a simplification of its elements. It was observed that the use of the complete Hessian matrix leads to different results; consequently, the simplified approach may not be recommended for some nonlinear parameter estimation problems, such as the one reported in this paper. Joint confidence regions played a key role in the analysis of parameter confidence intervals, especially in the order of the fractional differential equation. It is also important to conclude that the fractional model considered was in fact fractional, as the estimated order of the derivative was higher than its error and also statistically different from 1, showing that the fractional nature of the model is valid.



---

*Source: 642101-2013-05-21.xml*