Langmuir 1997, 13, 1303-1306
1303
Comparison of Various Numerical Procedures for Analysis of Structural Heterogeneity† A. M. Puziy,*,‡ V. V. Volkov,‡ O. I. Poznayeva,‡ V. I. Bogillo,§ and V. P. Shkilev§ Institute for Sorption and Problems of Endoecology, Kiev, Ukraine, and Institute for Surface Chemistry, Kiev, Ukraine Received December 20, 1995. In Final Form: December 24, 1996X Several methods (condensation approximation, regularization, and non-negative least-squares (NNLS)) have been applied for calculating the pore size distribution (PSD) function from simulated isotherms. It has been shown that two near by peaks in the low size region may be recovered well by both regularization and NNLS method only for isotherms without experimental error. Wide peaks may be recovered by numerical methods almost independently on a number of experimental points and error inherent in input data. Condensation approximation is the most insensitive to the experimental errors. However it gives a distorted one peak distribution shifted to low values of micropore size and having a long tail at large values. In the case of isotherms with large errors (10% for regularization and 5% for NNLS), both numerical methods give PSD functions close to that obtained by the CA method.
Introduction Adsorption on heterogeneous adsorbents plays an important role in many fields of science and technology. A quantitative characteristic of solid adsorbent with respect to its heterogeneity is the distribution function, which may be obtained using a variety of methods. For this purpose the adsorption measurements of standard substances from the gas phase are often used. Many methods have been developed for calculating the distribution function from adsorption data.1-3 Almost all of them deal with heterogeneity in terms of energetic distribution function, while for porous solids the structural heterogeneity, i.e., pore size distribution (PSD) function, is very important. In this paper we focus on the analysis of structural heterogeneity using several adsorption methods. Because there is no universal method for the characterization of solid adsorbents, we also wish to show a comparison of selected numerical methods for evaluation of structural heterogeneity using computer-simulated isotherms. In this study the following methods have been used: (i) condensation approximation (CA), very popular due to its simplicity and applied here as base reference method; (ii) regularization method; and (iii) non-negative least-squares (NNLS) methods. To enable a rigorous comparison, these methods were modified for analysis of the PSD function from adsorption data. Simulated Isotherms Since the theoretical description of physical adsorption on heterogeneous adsorbents is most often based on the adsorption integral equation, we used it here for simulation of benzene adsorption isotherms at 20 °C * Corresponding author: e-mail,
[email protected]; fax, (38-044) 293-06-39. † Presented at the Second International Symposium on Effects of Surface Heterogeneity in Adsorption and Catalysis on Solids, held in Poland/Slovakia, September 4-10, 1995. ‡ Institute for Sorption and Problems of Endoecology. § Institute for Surface Chemistry. X Abstract published in Advance ACS Abstracts, February 15, 1997. (1) Jaroniec, M.; Madey, R. Physical Adsorption on Heterogeneous Solids; Elsevier: Amsterdam, 1988. (2) Rudzinski, W.; Everett, D. H. Adsorption of Gases on Heterogeneous Surfaces; Academic Press: New York, 1992. (3) Cerofolini, G. F.; Re, N. Riv. Nuovo Cimento Soc. Ital. Fis. 1993, 16, 1.
S0743-7463(95)01560-5 CCC: $14.00
Θt )
∫∆Θl(P/P0,x) f(x) dx
(1)
where Θt is the normalized total adsorption isotherm as a function of relative pressure P/P0 at constant temperature, Θl is the “local” adsorption isotherm in pores with equal size x, f(x) is the PSD function, and ∆ is the range of pore size values (xmin ) 0.25 nm and xmax ) 2 nm in the case of benzene). As the local isotherm, the DubininAstakhov (DA) equation has been used with exponent n )3
[ ( )]
Θl ) exp -
A βE0
3
(2)
where A ) RT ln(P0/P), work of adsorption, β is the affinity constant, and E0 is the characteristic energy of adsorption, which is related to the halfwidth of micropores x4
E0 ) βk/x
(3)
The choice of parameter n ) 3 is based on the fact that truly homogeneous systems tend to follow the DA equation with a value of n close to 3.5-7 As the PSD function f(x) in eq 1, the tri-Gaussian distribution was used
f(x) )
∑i wi
1 σix2π
[
exp -
]
(x - x0)2 2σi2
(4)
with parameters listed in Table 1. From 10 to 100 points equally spaced along A ) RT ln(P0/P) were generated for a different P/P0 range and then random error was imposed on the adsorption value (Table 2). There was no error imposed on the P/P0 value. By simulating the above isotherms, we meant that the adsorption system benzene-active carbon is considered. The PSD function is represented by three peaks, which are characteristic for typical active carbons. Two sharp peaks at 0.5 and 0.7 nm, situated at the micropore region, were chosen to verify the ability of numerical procedures (4) Dubinin, M. M. Carbon 1985, 23, 373. (5) Stoeckli, H. F.; et al. Carbon 1989, 27, 125. (6) Stoeckli, H. F. Carbon 1989, 27, 962. (7) Stoecklii, H. F.; et al. Carbon 1990, 28, 907.
© 1997 American Chemical Society
1304 Langmuir, Vol. 13, No. 5, 1997
Puziy et al.
Table 1. Parameters of tri-Gaussian Distribution no.
weight of peak wi
mean x value (nm)
width of distribution, σi (nm)
1 2 3
0.05 0.05 0.90
0.5 0.7 1.1
0.02 0.02 0.25
Table 2. Simulated Isotherms no.
no. of points
P/P0 range
error, %
1 2 3 4 5 6 7
10-100 10-100 10-100 10-100 10-100 10-100 10-100
10-9-1.0 10-6-0.3 0.001-0.3 10-6-0.3 10-6-0.3 10-6-0.3 10-6-0.3
0 0 0 0.1 1 5 10
Heterogeneity Analysis Condensation Approximation. The most simple approximate method to obtain a distribution function from integral eq 1 is the condensation approximation (CA).2,8,9 Since this method is widely used for obtaining energy distribution and, to our knowledge, has not yet been applied for analysis of structural heterogeneity, we briefly show the general idea of the CA method in terms of the PSD function. According to CA the local isotherm is replaced by a onestep function
ΘCA ) 0 for A > Ec (5)
where A ) RT ln(P0/P) is work of adsorption and Ec is the critical characteristic energy of adsorption for micropores which are occupied at a given pressure. It is equivalent to an assumption that at a given relative pressure P/P0 those pores are filled in only, for which the critical energy of adsorption exceeds work of adsorption A. The step position is determined at Θ ) 0.5 by minimizing the square of the deviation between the true isotherm and step function. The relationship between halfwidth of micropores x and relative pressure P/P0, at which these pores are filled may be obtained by combination of the condensation at Θ ) 0.5 with DA eq 2 and equation 3
βk(ln 0.5)1/3 x)RT ln(P0/P)
(6)
Due to the fact that θCA is a step function, eq 1 may be simplified to
Θt )
∫∆f(x) dx
(7)
Then the PSD function may be obtained by differentiating eq 7
f(x) ) -dΘ/dx
∫∆Θlf(x) dx - Θt}2 + R∫∆f2(x) dx
Φ){
to recover two near by peaks. A third wide peak is characteristic for highly activated carbons.
ΘCA ) 1 for A e Ec
Regularization Method. From a mathematical point of view, eq 1 is a Fredholm integral equation of the first kind. Solving this equation with respect to the PSD function is a numerically ill-posed problem; i.e., small changes in measured adsorption Θt caused by experimental errors can distort significantly the sought for function. One of methods to solve ill-posed problems is the regularization method.10,11 The fundamental idea of the numerical regularization method is to replace the illposed problem of minimizing the select function by a wellposed problem which smoothes the calculated distribution function and distorts the origin problem insignificantly. Hence, solving the integral equation is replaced by minimization of the functional
(8)
Thus, the experimental adsorption isotherm Θt, measured as function of relative pressure P/P0, may be represented as a function of pore half-width x using eq 6 and then differentiated to obtain the desired PSD function. (8) Roginskij, S. S. Dokl. Akad. Nauk USSR, Ser. B 1944, 45, 194. (9) Harris, L. B. Surf. Sci. 1968, 10, 129.
(9)
where R is a regularization parameter (0 < R e 1). This method was successfully employed for obtaining the energy distribution function from gas adsorption data.12-16 However, it has not been yet widely used for analysis of structural heterogeneity, that is, for finding the PSD function.17,18 The choice of the optimal of regularization parameter R is crucial for the sought for PSD function. A very low value of regularization parameter gives rise to spurious peaks; while too high a value oversmooths the PSD function. A detailed description of strategies for finding the optimal R value in adsorption application is given in refs 12 and 19. In the present paper the following method, consisting of two steps, has been used for finding the optimal regularization parameter. First, the functional (9) is minimized at R ) 0. Then the value of this minimum may serve as measure of the accuracy of experimental data
[∑
]
1 m (Θt - Θlf0(x))2 r ) min f mi)1
1/2
(10)
In the second step, the value of R is calculated at which the PSD function f1(x) minimizing the functional (9) satisfies the following condition:
[∑ 1
m
mi)1
]
1/2
(Θlf0(x) - Θlf1(x))2
) βr
(11)
In the present work parameter β was taken as equal to 1. The special convenience of this method is that the solution may be improved by varying the free parameter β. The necessity of the free parameter stems from the fact that there does not exist a method of finding the stable solution for ill-posed problems without taking into account errors in input data.20 This means that if accuracy of (10) Tikhonov, A. N. Dokl. Akad. Nauk USSR 1943, 39, 195. (11) House, W. A. J. Colloid Interface Sci. 1978, 67, 166. (12) Szombathely, M.v.; Brauer, P.; Jaroniec, M. J. Comput. Chem. 1992, 13, 17. (13) Heuchel, M.; Jaroniec, M.; Gilpin, R. K. Langmuir 1993, 9, 2337. (14) Heuchel, M.; Jaroniec, M. Langmuir 1995, 11, 1297. (15) Heuchel, M.; Jaroniec, M. Langmuir 1995, 11, 4532. (16) Papenhuijzen, J.; Koopal, L. K. In Adsorption from Solutions; Ottewill, R. H., Rochester, C. H., Smoth, A. L., Eds.; Academic Press: London, 1983; p 211. (17) Puziy, A. M.; Leboda, R.; Bogillo, V. I.; Shkilev, V. P.; Lodyga, A. Adsorpt. Sci. Technol. 1995, 12, 267. (18) Leboda, R.; Skubiszewska-Zieba, J. S.; Bogillo, V. I. Langmuir 1997, 13, 0000. (19) Hansen, C. SIAM Rev. 1992, 34, 561. (20) Leonov, A. S.; Jagola, A. G. Vestn. Mosk. Univ., Ser. 3: Fiz., Astron. 1995, 36, 28.
Analysis of Structural Heterogeneity
Langmuir, Vol. 13, No. 5, 1997 1305
measured isotherm is unknown, the choice of a regularization parameter will always be subjective. None of the existing methods for choosing parameter R can be strictly proved and all of them give, in one case or another, unsatisfactory results. Without a free parameter β one has to be satisfied with the solution even if it is either oversmoothed or contains spurious peaks. The applied computer program was used here for calculation of the PSD function with respect to positive values. In this program the integrals were computed by the trapezoid method. Minimization of the functional (9) was performed by the conjugate gradient method.21 The advantage of the regularization method is the possibility to calculate more points on a desired PSD function than the number of experimental points available. Non-negative Least-Squares Method. A second numerical method used here to solve the integral equation (1) is the non-negative solution of a system of linear equations (NNLS method).22 In this method the integral equation (1) is represented by set of linear equations:
Θt ) Af
(12)
where Θt is measured adsorption, f is the PSD function f(x), and A is a matrix with elements ai,j ) Θl (P/P0,x). The original NNLS algorithm implies the search for the largest possible set of columns from the original matrix A which still gives only positive values for the unknown component of vector f by the consecutive iterative procedure. Thus the NNLS algorithm calculates the trial matrix Ap which as well as original matrix A is also near singular. To overcome the instability in inverting near singular matrix Ap the special procedure has been included in the original NNLS algorithm. This procedure includes reduction of matrix Ap to upper triangular form and replacing rows having a diagonal value lower than regularization parameter τ by 0.22 The larger the value of τ used, the more regularized the solution may be obtained. The regularization parameter τ has been chosen to satisfy the following condition
[
1
∑ mi)1
]
1/2
m
(Θt - Θt,c)2
)σ
(13)
where Θt and Θt,c are the experimental and calculated isotherms, m is the number of experimental points on the isotherm, and σ is the error inherent in experimental data. The advantages of the NNLS procedure are its stability and finiteness. The shortcomings are the relatively large number of iterations and presence of spurious zero elements inside the interval of positive values. This creates oscillating behavior of the PSD function, which obviously has no physical meaning. Also a limitation of this method is that the number of points on the desired PSD function cannot exceed the number of experimental points involved in the calculation. In the present paper the following improvement has been proposed to the standard NNLS algorithm. The standard NNLS procedure implies calculation of the dual matrix, which serves as a criterion for introducing a new column. Instead of using this criterion (optimal choice), we used a random choice. Thus, dual matrix needs not be calculated on each iteration. This decreases the computational time and allows to obtain the same solution for a lesser number of iterations (Table 3). Also, random choice eliminates spurious zero elements in the PSD (21) Bogillo, V. I.; Shkilev, V. P. Langmuir 1996, 12, 109. (22) Lawson, C. L.; Hanson, R. J. Solving Least Squares Problems; Prentice-Hall: Englewood Cliffs, NJ, 1974.
Figure 1. PSD functions obtained by the CA method from simulated isotherms abbreviated as: number of isotherms from Table 2; number of points. Table 3. Comparison of Random and Optimal Choice in NNLS Method for 50-Point Isotherm no. of iterations no. of isotherm
random
optimal
1 3 5 6
140 86 103 98
328 245 330 153
function. The optimal choice is preferable only in cases of small content of nonzero elements. Results and Discussion Since the CA method uses a first derivative isotherm, problem with spurious peaks may arise, if noisy data will be differentiated without preprocessing. Examples of preprocessing may be a smoothing spline technique23 or using a DA isotherm as smoothing equation.24 Pore size distribution functions obtained by the CA method are presented in Figure 1. As pointed out earlier,2 the CA method gives a flattened peak. Small sharp peaks are not recovered at all, even if an isotherm free of experimental error is used (Figure 1, panel a). The maximum of the resulting PSD function is shifted to lesser values of halfwidth, obviously due to the presence of two small peaks on original distribution. If a proper smoothing routine is applied, the CA method is insensitive to experimental error. Even as large as 10% error (Figure 1, panel c) has a very small impact on resulting PSD function as compared with that obtained from an isotherm without experimental error (Figure 1, panel a). The quantity of experimental points is of no significance for the CA method: a 10-point isotherm gives the same PSD as a 100-point isotherm (Figure 1, panels b and c). The CA method gives in fact quite the same PSD function from all isotherms. Due to the fact that in the CA method the isotherm is approximated by a step function, the range of PSD depends (23) Nederlof, M. M.; et al. Environ. Sci. Technol. 1994, 28, 1037. (24) Puziy, A. M. Langmuir 1995, 11, 543.
1306 Langmuir, Vol. 13, No. 5, 1997
Puziy et al.
Figure 3. PSD functions obtained by the NNLS method. The explanation of curves is given in Figure 1.
Figure 2. PSD functions obtained by the regularization method. Explanation of curves as in Figure 1. Table 4. Halfwidth of Micropores Occupied at Given Relative Pressure (Calculated for Benzene at 293 K) P/P0
x, nm
P/P0
x, nm
2.67 × 10-8 10-6 10-5 10-4
0.25 0.32 0.38 0.47
0.001 0.01 0.100 0.113
0.63 0.95 1.89 2.00
on the interval of measured pressure. To calculate the PSD function in the entire range of micropores (x ) 0.25-2 nm), measurements of the adsorption of benzene in the range of P/P0 ) 2.67 × 10-8 to 0.113 should be available. This range may easily be calculated by eq 6 (Table 4). Contrary to the CA method, there is no such limitation when numeric methods are used. The lowest limit of micropore size in this case may be obtained even if only a “window” of data is available. Figure 2 and Figure 3 present PSD functions obtained by the regularization and NNLS methods. Both numerical methods give similar results. Sharp peaks are recovered well only for isotherms with a wide range of relative pressure and without experimental errors (isotherms 1 and 2, see Table 2). Eliminating from calculation the points measured at low relative pressure (isotherm 3 simulated for P/P0 ) 0.0010.3) as well as low experimental error (isotherms 4 and 5 simulated with 0.1 and 1% error) gives only a shoulder instead of peaks at 0.5 and 0.7 nm.
Although a wide peak at 1.1 nm is recovered quite well by both numerical methods, the regularization is more insensitive for large errors. Only a 10% error on isotherm can distort the PSD function obtained by regularization (Figure 2, panel d), while for PSD obtained by the NNLS method, this distortion is observed already at 5% error (Figure 3, panel c). Isotherms with large experimental error give a PSD function close to that obtained by the CA method. Regularization recovers well the PSD function independently on number of experimental points on the isotherm (Figure 2, panel a), while for the NNLS method it is not possible to obtain good results for 10-point isotherms. Conclusions Two near by peaks on the low size region may be recovered well by both regularization and NNLS methods only for isotherms without experimental error. Even error as low as 0.1% flattens sharp peaks and gives only a shoulder. Wide peaks may be recovered by numerical methods almost independently on number of experimental points and error inherent in input data. Condensation approximation is the most insensitive to the experimental errors. However it gives a distorted one peak PSD shifted to low values of micropore size with a long tail at large values. In the case of isotherm with large errors (10% for regularization and 5% for NNLS), both numerical methods give a PSD function close to that obtained by the CA method. LA951560S
