Solution of an inverse adsorption problem with an epidemic genetic algorithm and the generalized extremal optimization algorithm

Abstract

In the present work two recently developed stochastic methods, the epidemic genetic algorithm and the generalized extremal optimization algorithm are used for the solution of an inverse mass transfer problem, which is implicitly formulated as an optimization problem, for the estimation of parameters associated with the adsorption of biomolecules in resin beds. The estimates obtained with both methods present good accuracy, even in the presence of noisy data, provided that the model and experiment used are sensitive to the parameters being estimated. With Thomas' model for the direct mass transfer problem and real experimental data for lisozyme in adsorption columns, it is possible to estimate the maximum adsorption capacity in Langmuir's adsorption isotherm.

Keywords:

• adsorption
• optimization
• Thomas' model
• inverse problems

1. Introduction

Due to the relevant applications in the food and pharmaceutical industries, there is a growing demand for the formulation and solution of direct and inverse mass transfer problems 1–5. A high demand for purification processes for increasingly complex substances has been observed in recent years, with one of the most promising alternatives being the simulated moving bed (SMB) chromatography. For a full understanding of the operation and optimization 6 of SMBs, and a possible scale-up to industrial production, an accurate knowledge of mass transfer mechanisms and their dependence on the physico-chemical and process parameters involved is required. The first step in that direction consists of the characterization of adsorption columns 7.

Kowalczyk et al. 8 estimated the pore-size distribution function that is used as a quantitative characteristic of the porous structure of solid adsorbents with respect to their heterogeneity, and Kowalczyk et al. 9 used a simple adsorption genetic algorithm (GA) for the estimation of parameters in a equation that models the adsorption process in a homogeneous micropore system.

The estimation of adsorption parameters for chromatography systems is a very important step in the column characterization used in the design of continuous separation processes based on the adsorption phenomenon, such as in simulated moving beds. For that purpose Forssén et al. 10 and Câmara et al. 11 used the least squares approach, the former having calculated the Jacobian matrix with numerical differentiation with complex variables. Vasconcellos et al. 12 used GAs. Lage et al. 13 have also used GA and a hybridization of GA with the Levenberg-Marquardt method (GA-LM), while Lage et al. 14 applied artificial neural networks (ANNs) developed for the solution of the inverse mass transfer problem, and Lage et al. 15 used a hybridization of ANNs with the Levenberg-Marquardt method (LM).

Yu et al. 16 and Ziyang et al. 17 have estimated adsorption equilibrium constants and kinetic parameters using GA. James et al. 18 estimated adsorption isotherms coefficients with the conjugate gradient method, and Zhang et al. 19 used a variation of GA for the same purpose.

In the present work, two stochastic methods are applied for the estimation of adsorption isotherm coefficients: (i) the generalized extremal optimization (GEO) 20; and (ii) the epidemic genetic algorithm (EGA) 21.

For the solution of the direct problem Thomas' model is used, which provides accurate results when the effects of axial dispersion are negligible in comparison with other mass transfer mechanisms 7, 2223.

2. Mathematical formulation and solution of the direct problem

A mobile liquid phase composed by a diluted solution of the adsorbate of interest, for example a biomolecule, percolates through a resin bed, the solid fixed phase. The adsorbate is transferred from the bulk of the solution to the vicinity of the resin particles, i.e. a mass transfer mechanism through the liquid film, and then it diffuses to the interior of the particle pores being then adsorbed by the solid matrix.

The mathematical model for the separation chromatographic process is based on the mass balance for the two phases; one for the mobile phase that flows through the macroscale porous fixed resin bed, and the other for the resin particles involving the microscale porous solid matrix 24.

For the particular case in which the axial dispersion may be neglected, and the adsorbate inflow to the column of resins is constant, with concentration C₀, the adsorption problem described here has an analytical solution 25–27, which was first derived by Thomas 28. According to this solution, the adsorbate concentration in the mobile liquid phase at the exit of the adsorption column, C, as a function of time, t, i.e. the breakthrough curve, is given by (1) where (2) (3) (4) Q is the volumetric flow rate, h is the total length of the adsorption column, A_c is the cross section of the column, k_d is the dissociation rate constant defined as (5) k₁ is the adsorption rate constant, k₂ is the desorption rate constant and q_m is the maximum adsorbate concentration that the adsorbent may adsorb (maximum adsorption capacity), according to the Langmuir adsorption isotherm, (6) The function J in Equation (1) corresponds to (7) where I₀ is the modified Bessel function of the first kind and order zero, which may be approximated by an asymptotic series whose two first terms are given by (8) When the geometry, initial condition, boundary conditions, physico-chemical and process parameters are known, we are able to calculate the concentration of the adsorbate at the exit of the column, i.e. C(h,t). This is the direct problem, i.e. the determination of the breakthrough curve.

3. Mathematical formulation and solution of the inverse problem

Here, we are interested in the estimation of the coefficients in the adsorption isotherm given by Equation (6) using measured values of the adsorbate concentration at the exit of the column, i.e. C_measi = C_meas(t_i) with i = 1,2,…,N_d, where N_d represents the total number of experimental data available.

According to the sensitivity analysis performed by Folly et al. 7 we may try to estimate the vector of unknowns (9) but difficulties may arise due to low values of the sensitivity coefficients related to the parameter k_d, and that behaviour was in fact observed by Lage et al. 15.

In the present work besides the estimation of we will also look into the estimation of (10) fixing the value of k_d, and also into the estimation of (11) fixing the value of q_m.

As the number of experimental data available, N_d, will be much larger than the number of unknowns considered, N_u = 1 or 2, we formulate the inverse problem implicitly as an optimization problem in which we seek to minimize the cost function given by the summation of the squared residues between the calculated and measured adsorbate concentrations at the exit of the adsorption column, i.e. C_calci and C_measi, respectively, with i = 1,2,…,N_d, (12) In order to find the value of for which S is minimum, we have used two recently proposed stochastic methods: (i) GEO 20; and (ii) EGA 21. Both the methods will be described next.

4. The generalized extremal optimization algorithm

GEO has been proposed and implemented recently 20, and was devised to be applied in complex optimization problems. Based on the Back-Sneppen simplified model of evolution 29, it has been applied successfully to design optimization problems 30, 31 as well as to an inverse radiative transfer problem 32.

GEO makes no use of derivatives requiring only the solution of the direct problem given by Equations (1–4). The main steps of the canonical GEO algorithm are shown in Figure 1.

Figure 1. The canonical GEO algorithm.

A very attractive feature of GEO is that it has only one free parameter to adjust, τ, which may be considered an advantage when compared to other stochastic methods that have a larger number of parameters to be set.

In GEO, a string of L bits, which encodes the N_u unknowns, is considered a population of species, i.e. each bit represents a species. Each bit is associated to a fitness number that is proportional to the gain, or loss, the cost function value has in flipping that particular bit. All bits are then ranked from k = 1, for the least adapted bit, to k = L for the best adapted. A bit is then mutated according to the probability distribution , where k is the rank of a selected bit candidate to mutate. Making all bits have the same probability to mutate, whereas for only the worst adapted bit will mutate.

In practice, it has been observed that the best value of τ, i.e. the one that yields the best performance of the algorithm for a given application, usually lies in the range (0.75,3.0). A detailed description of GEO can be found in Sousa et al. 20.

5. The epidemic genetic algorithm

Simple GAs operate on a fixed-sized population of fixed-length individuals, and in general the individuals are represented by a binary string that encodes the variables of the problem that the algorithm is trying to minimize.

In the present work, each individual is composed by the unknowns of the inverse problem, which are given by Equations (9–11), and in fact is represented by a real-valued string.

Simple GAs use basically three operators: selection, crossover and mutation. The selection operator identifies the fittest individuals of a given population to serve as parents of the next generation. The fitness value of each individual is related to the value of the cost function given by Equation (12). The selection operator ensures that the best fit individuals have a higher probability to be selected to reproduce and form a new generation.

The crossover operator randomly chooses a pair of individuals among those previously selected to breed and exchanges information between them.

The mutation operator is usually considered a secondary operator. Its main function is to restore diversity that may be lost from the repeated application of the selection and crossover operators. This operator simply takes one string from the population and randomly alters some value within it. Following the example of nature, the probability of applying the mutation operator is very low compared to the probability of applying the crossover operator.

In the present work, a fourth operator is applied, the epidemical strategy 21, leading to EGA. This operator is activated when a prescribed number of generations are reached without improvement of the best individual. Then, all the population is affected by a plague, and only those that have the best fit, say the top 10% individuals better fitted, survive. The remaining individuals die and are replaced by new individuals with new genetic variability, such as immigrants arriving in order to evolve the population. Two parameters need to be chosen: one determines when the strategy will be activated, i.e. the number of generations without improvement of the fitness value of the best individual, N_e, and the other determines the fraction of the population that will survive the plague, f_e.

Slow convergence is a particular problem for the stochastic optimization schemes, the same observation is valid to the GA. Alternative GA techniques have been designed to improve the convergence performance. Another challenge is to design a GA for dealing with multimodal function optimization. For attending to such features some genetic operators were proposed such as: micro-GA 33 and niching 34. A niche is a particular sub-domain of the entire search space, where each niche has different nuances on each other. This is a natural feature to maintain the genotype diversity in a population 35–37, under insulation these niches are a mechanism for the speciation. Goldberg and Richardson 34 have introduced the method of sharing functions–now the method is known as sharing, and it is directly applicable to the multimodal function optimization 36. Sharing is a natural algorithm to implement niching 38.

Another procedure to preserve the diversity in the GA population is the re-initialization. Mahfoud 36 has reported that Goldberg 39 has investigated such ideas, but re-initialization is not treated as a new GA operator. Micro-genetic algorithms are GAs with small populations with re-initialization. Krishnakumar's 33 micro-GA approach uses a population with five individuals, driven by elitism. He compared his micro-GA with standard GA (SGA), resulting in a faster and better solution to a engineering control problem. Other applications of micro-GA have been presented by Liu and Han 37. We pointed out that Chakraborti et al. 40 addressed a solution, an important and difficult real-world problem, where it was shown to be of little additional advantage over the simple GA, when the objective function is stationary. Micro-GA deals with selection and crossover GA operators, while mutation is usually omitted 36, 37. On the other hand, EGA employs all standard GA operators, and also the typical population size of the SGA, but it introduces the epidemic strategy.

An interesting issue would be a comparison among the diversity preservation strategies (niching, micro-GA, EGA, for example). This will be the focus of a future work. Here, the goal is just the performance comparison between GEO and EGA.

6. Results and discussion

Because of the availability of real experimental data, we have chosen the system investigated by Chase 26 for the substance lisozyme as our test case. Table 1 presents the process and physico-chemical parameters used.

Table 1. Process and physico-chemical parameters for Chase's experiment with lisozyme 26.

Parameter	Value
H–column height (cm)	10.4
Ac–column cross-section (cm^2)	0.785
Q–volumetric flow rate (mL min^−1)	1.0
C0–adsorbate concentration at the column inlet (mg mL^−1)	0.1
qm–maximum adsorption capacity (mg mL^−1)	14
kd–dissociation rate constant (mg mL^−1)	0.025
k1–adsorption rate constant (mL mg · min^−1)	0.20

The values of k₁, k_d and q_m shown in Table 1 were obtained using a batch experiment and considering a Langmuir adsorption isotherm.

From the sensitivity analysis performed by Folly et al. 7 one concludes that using Thomas' model and the experimental data obtained by Chase 26 for lisozyme it is not possible to estimate k₁. Therefore, in all computations, the results of which will be presented next, a fixed value for this parameter has been considered, i.e. the one shown in Table 1.

Tables 2 and 3 show the estimates obtained for the vector of unknowns , see Equation (9), for lisozyme, using GEO and EGA, respectively, and Chase's experimental data 26. Each table shows the results obtained for 10 different runs of each algorithm, as well as the average and the standard deviation for each unknown. It also shows the value of a measure of the dispersion of the estimates given by (13)

Table 2. Estimates for k_d and q_m for lisozyme using GEO with τ = 1.25 and Chase's experimental data 26. d = 0.11514 (Equation (13)).

Run	kd (mg mL^−1)	qm (mg mL^−1)	S Equation (12) (mg mL^−1)^2
1	0.0272	11.2862	0.0251
2	0.0335	11.8862	0.0243
3	0.0350	14.3182	0.0233
4	0.0287	14.0856	0.0327
5	0.0258	13.2334	0.0278
6	0.0325	12.0309	0.0231
7	0.0321	12.7347	0.0254
8	0.0215	12.7701	0.0236
9	0.0325	13.9952	0.0234
10	0.0316	11.0701	0.0235
Average μ	0.0300	12.7411
Standard deviation σ	0.0042	1.1669
	13.82	9.16

Table 3. Estimates for k_d and q_m for lisozyme using EGA with an initial population of 100 individuals, Chase's experimental data 26 and n_e = 3 and f_e = 0.1. d = 0.047396 (Equation (13)).

Run	kd (mg mL^−1)	qm (mg mL^−1)	S Equation (12) (mg mL^−1)^2
1	0.0372	14.0389	0.0240
2	0.0455	14.1111	0.0369
3	0.0475	14.2639	0.0270
4	0.0448	14.1780	0.0210
5	0.0446	14.1892	0.0205
6	0.0444	14.1681	0.0205
7	0.0406	14.0509	0.0208
8	0.0397	14.0909	0.0221
9	0.0374	13.9978	0.0241
10	0.0395	14.0770	0.0215
Average μ	0.0421	14.1166
Standard deviation σ	0.0037	0.0817
	8.74	0.58

Figures 2 and 3 show graphically the dispersion of the results obtained with the stochastic methods GEO and EGA, respectively. Observe that the scales used in the two figures are different.

Figure 2. Dispersion of the results obtained with GEO.

Figure 3. Dispersion of the results obtained with EGA.

Figure 4 shows Chase's experimental data 26 and the breakthrough curves calculated using the average values obtained for q_m and k_d, shown in Tables 2 and 3, using GEO and EGA, respectively.

Figure 4. Chase's experimental data 26 and calculated breakthrough curves using the average of the estimated values for q_m and k_d with GEO and EGA (Tables 2 and 3).

From the results presented we conclude that the low sensitivity to the parameter k_d may be affecting negatively the estimation of q_m when the two parameters are estimated simultaneously.

We may also conclude that with the control parameters chosen for GEO and EGA the latter presents a lower dispersion of the estimates. This fact is confirmed by the value of the measure d, given by Equation (13), shown in Tables 2 and 3.

We have then fixed the value of q_m as 14 mg mL⁻¹ and have estimated k_d using both GEO and EGA. The results obtained after 10 runs of both algorithms are shown in Table 4.

Table 4. Estimated values for k_d, being q_m fixed at 14 mg mL⁻¹, after 10 runs of GEO and EGA.

	GEO	EGA
(mg mL^−1)	0.0261	0.0381
σkd (mg mL^−1)	0.0009	4.438×10^−5
	3.34	0.12

Even though the dispersion of the estimates is smaller with EGA we observe from Tables 1 and 4 that the estimated average value for k_d obtained with GEO is closer to the experimental value shown in Table 1. Nonetheless, as mentioned before, due to the low value of the sensitivity coefficient related to this parameter the estimates may not be accurate.

Figure 5 shows Chase's experimental data and the breakthrough curves calculated using the values for q_m and k_d shown in Table 4.

Figure 5. Chase's experimental data 26 and calculated breakthrough curves using the average estimated values for k_d with GEO and EGA being q_m fixed at 14 mg mL⁻¹ (Table 4).

Finally, we have varied the value of k_d in the range from 0.015 mg mL⁻¹ up to 0.035 mg mL⁻¹, and for each fixed value of k_d we have estimated q_m using both GEO and EGA. The results obtained after 10 runs of both algorithms are shown in Tables 5 and 6.

Table 5. Estimated values for q_m, being k_d fixed, after 10 runs of GEO.

kd (mg mL^−1)	(mg mL^−1)	σqm (mg mL^−1)
0.015	14.92	0.148	0.99
0.020	14.52	0.137	0.94
0.025	14.10	0.184	1.31
0.030	14.01	0.153	1.09
0.035	13.99	0.139	1.00

Table 6. Estimated values for q_m, being k_d fixed, after 10 runs of EGA.

kd (mg mL^−1)	(mg mL^−1)	σqm (mg mL^−1)
0.015	15.0	0	0
0.020	14.507	0.0035	0.02
0.025	14.201	0.0014	0.01
0.030	14.046	0.0006	0.005
0.035	14.007	0.0004	0.003

Figures 6 and 7 show Chase's experimental data 26 and the calculated breakthrough curves for lisozyme with the fixed value of k_d and the estimated values of q_m with GEO and EGA, respectively.

Figure 6. Chase's experimental data 26 and calculated breakthrough curves using the fixed value of k_d and the average estimated values for q_m with GEO (Table 5).

Figure 7. Chase's experimental data 26 and calculated breakthrough curves using the fixed value of k_d and the average estimated values for q_m with EGA (Table 6).

From the results shown in Tables 5 and 6 we observe a good agreement with the experimental value of q_m shown in Table 1 and again it seems that the dispersion of the estimates is smaller with EGA.

From Figures 6 and 7 it is observed that k_d, as expected, seems to have a negligible effect on the breakthrough curves.

All results shown for GEO were obtained using 5000 evaluations of the cost function. The dispersion of the estimates may become smaller if a larger number of function evaluations is considered. Of course a higher amount of CPU time will be required. This subject must be further investigated.

Due to the simplicity of the solution obtained for the direct problem with Thomas' model, only 10 generations were required to obtain the results presented before using EGA. Tests were performed with a total of 100 or 1000 generations, nonetheless no improvement was observed in the value of the cost function.

Concerning the accuracy of the solutions found by GEO and EGA when compared to the experimental data, for the results with fixed k_d, both algorithms presented similar average results, with the EGA ones presenting a lower dispersion. On the other hand, when estimating k_d and q_m simultaneously, the EGA had clearly a better performance than GEO.

7. Conclusions

The application of two stochastic methods, EGA and GEO yielded good estimates for the maximum adsorption capacity, one parameter in Langmuir's adsorption isotherm, for a system with lisozyme in a diluted solution being adsorbed in resin beds.

It was observed that EGA tends to provide estimates with a smaller dispersion when compared to GEO. Nonetheless this conclusion cannot yet be generalized. A small number of function evaluations was considered for GEO. With a larger number of functions evaluations a reduction in the dispersion may occur.

The results obtained so far are very encouraging and the application of the recently developed stochastic methods in the solution of inverse mass transfer problems deserves further investigation.

Acknowledgements

The authors acknowledge the financial support provided by CNPq, Conselho Nacional de Desnvolvimento Científico e Tecnológico, and FAPERJ, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro.