A Metaheuristic Tabu Search Optimization Algorithm: Applications to Chemical and Environmental Processes

Abstract

Stochastic optimization methods are increasingly used for optimizing processes that are difficult to solve by conventional techniques. These methods are widely employed to optimize the processes which have higher dimensionality with severe nonlinearities. Different methods of this kind include the genetic algorithm (GA), simulated annealing (SA), differential evolution (DE), ant colony optimization (ACO), tabu search (TS), particle swarm optimization (PSO), artificial bee colony (ABC) algorithm, and cuckoo search (CS) algorithm. Among these methods, tabu search (TS) is a potential tool used to find a feasible optimal solution from a finite set of solutions. The memory used in TS will remember the current best solution and it also enables the TS to track the last solutions while guiding the search moves. The capability of memory and strategic adaptation features of TS enable it to make use of good solutions and also search for new feasible regions in the search space. TS has been successfully applied to solve a wide spectrum of optimization problems in different disciplines. This chapter describes the TS algorithm in detail and its applications to chemical and environmental processes, specifically, dynamic optimization of a copolymerization reactor and inverse modeling of a biofilm reactor. In dynamic optimization of copolymerization reactor, the meta heuristic Tabu search (TS) is designed and applied to determine the optimal control policies of a styrene–acrylonitrile (SAN) copolymerization reactor. In inverse modeling of biofilm reactor, the tabu search is designed and applied to determine the parameters of kinetic and film thickness models as consequence of the validation of the mathematical models of the process with the aid of measured data acquired from an experimental fixed bed anaerobic biofilm reactor used in the treatment of pharmaceutical industry wastewater. For both the cases, optimization by Tabu search is carried out by suitably formulating the desired objective functions and the problems are solved by encoding the variables and parameters using real floating point numbers. The results explain the efficacy of TS for optimal control of polymerization reactor and inverse modeling of biofilm reactor.

1. Introduction

Tabu search is a meta heuristic problem solving approach used to solve combinatorial optimization problems. It was first proposed by Glover [1] and further developed by Hansen [2]. TS has now become an established search procedure and has been successfully applied to solve a wide spectrum of optimization problems [3, 4, 5, 6, 7, 8, 9].

1.1 Basic principle

TS uses a memory which allows it to remember the current best solution and also enables it to explore the previous solutions and to direct the search moves. The features of memory adaptation and exploration facilitate TS to find better solutions and discover new potential regions in the search space. The memory adaptation feature helps to realize its course of action to exploit the search space efficiently. The ability of TS to explore good solutions enables it to assimilate the intelligent search mechanisms to search for good solutions and to discover new potential regions. The use of adaptive memory helps TS to learn and creates a more flexible and effective search strategy compared with memory less methods, such as simulated annealing (SA) and genetic algorithms (GA).

1.2 Components of TS

The components of TS are explained as follows.

1.2.1 Neighbors generations and neighborhood search

To optimize the function f(x) globally from all the probable solutions x∈X in the space X, it requires to specify a structure in the vicinity of the solution space and the staring solution. The search advances to alter the existing solution to create a set of promising solutions in the vicinity of solution space. During the search process, the number of solutions traversed by TS is the product of the number of solutions in the vicinity of the solution space, N(x_i), and the number of iterations, k. The function at each iteration is evaluated for N(x_i) solutions. The better move is chosen in the vicinity of the complete solution space. The search proceeds to the next iteration, to find a solution in the vicinity of the accepted move. Thus, the TS builds a set of viable solutions by using a history record of the search.

1.2.2 Tabu list

The tabu list in TS maintains the data of the previously visited solutions. The list is included with the latest moves and it is altered dynamically as the search proceeds. The data in the tabu list helps to direct the move from the present solution to the next solution. At each iteration, the search process is maintained by updating the tabu list. The tabu list also avoids re-visiting of recent neighbors recorded in the list and thus save computational time.

1.2.3 Short-term memory and long-term memory

The information is stored in tabu list as recency-based short-term memory (RSM) and a frequency-based long-term memory (FRM). When the search proceeds, the nearly better solution to the present solution is sorted as tabu, and added to the recency-based tabu list. As fresh solutions enter the list, older solutions are removed from the bottom. Long-term memory relies on the frequency of a solution that is visited. When the tabu list is filled with the highest number of elements in the frequency based tabu list, the solution with the smallest frequency index will be replaced.

1.2.4 Intensification and diversification

These strategies are used to create neighbors that have higher likelihood of finding optimal solutions based on the data in the tabu list. Intensification strategies are used to search potential areas in detail around the areas that are found good in the past. Intensification strategies are generally employed based on long term memory whose components are used to create neighbors for search intensification [4, 6, 10]. Diversification strategies are used to search the complete viable region, thus restricting the search getting trapped in local optima. These strategies promote probing unvisited regions by creating solutions radically different from those searched earlier [4]. A frequency-based tabu list is used to keep track of the search area.

During the generation of neighbor solutions, the difference between the present and fresh neighbor solution is managed by using a coefficient, α. The change from the current point is multiplied by α during the course of building new neighbors. The coefficient α is in the form of a sine function [10].

$$\alpha =\frac{1}{2}\left[1+sin\left[\frac{\mathit{i\theta \pi }}{N_{\mathit{\text{neigh}}}}\right]\right]$$E1

Here i is index of the neighbor, N_neigh is the total number of neighbor solutions generated at each iteration, and θ is a parameter that controls the oscillation period of α.

1.2.5 Aspiration criterion

The TS conditions at times prevent moves leading to unvisited solutions. Aspiration criterion is a condition that can override the tabu status of a certain move. To avoid certain missing solutions during the search, the aspiration criterion in certain cases may invalidate the tabu property and maintains an appropriate balance between diversification and intensification [4, 10, 11].

An aspiration criterion that is designed based on a sigmoid function as given by

$$S\left(k\right)=\frac{1}{1+e^{-\sigma \left(k-k_{\mathit{\text{center}}}\times M\right)}}$$E2

where k_center, k, σ and M denote the tuning parameter, current iteration number, another tuning parameter and maximum number of iterations. The value of k_center can be in the range of 0.30–0.70, the value of σ can be in the range of 5/M to 10/M. A random number P that lies between 0 and 1 is generated from a uniform distribution at each iteration. If P is greater than S(k), the tabu property is active and the best non-tabu neighbor is used as a fresh starting point. If P is less than or equal to S(k), the aspiration criterion ignores the tabu property.

1.2.6 Stopping criteria

A stopping criterion is needed to terminate the search when the optimum is reached. The stopping criterion can be in the form of fixing the number of iterations or specifying a threshold for convergence of solution. The criteria like the maximum time termination [6] and termination-on convergence criteria [10] are also used as stopping criteria for the search process. The termination-on-convergence criteria is expressed by Lin and Miller [10] as.

$$\left|\frac{f_k\left(x\right)-f_{k-\mathrm{\Gamma }}\left(x\right)}{f_{k-\mathrm{\Gamma }}\left(x\right)}\right|<\delta$$E3

where δ is the ratio of the change in the objective function value, Γ = ηM, and η is the fraction of the maximum iterations (M) by which the change in the objective function is compared. As per this stopping criterion, if the enhancement over Γ generations is no longer than a threshold (δ), continuation of further iterations can be ineffective, and the search should be discontinued.

1.3 TS implementation procedure

The flow chart of TS algorithm is shown in Figure 1. Tabu Search begins with an initial solution x_o. The neighbor solutions are created by altering the existing solution through a sequence of moves. The best new neighbor, x^*, is used as the starting point for the next iteration, unless it is in the tabu list. Thus, even if no neighbor solutions are better than the starting solution, the best solution is still selected as the starting point for the next iteration. A record of the best solutions ever found, x^*, is separately maintained. Also, the adaptive memory in tabu lists guides the search by utilizing the benefit of past information. This memory facilitates TS to make strategic choices and accomplish responsive exploration.

Figure 1.

Flow chart of the TS algorithm.

TS is implemented using the following steps.

1. Starts with an initial solution x_o

2. Find the best solution, x^* and define tabu list

3. Neighbor’s generations and neighborhood search

   1. Define a neighborhood structure with respect to the solution space

   2. Change the present solution to form a neighborhood of possible solutions

   3. TS explore neighbor solution at each iteration

   4. Searches entire neighborhood and select enhanced solution

4. Tabu list

   1. Majority of the recent solutions are in tabu list

   2. Older solutions are discarded

5. Short-term memory and long-term memory

   1. Recency-based short-term memory (RSM)

   2. Frequency-based long-term memory (FRM)

6. Intensification and Diversification

   1. To search potential areas

   2. To search whole feasible region

7. Aspiration criterion

   1. To explore unvisited solutions

   2. To evade feasible missing solutions

8. Stopping criteria

   1. According iterations

   2. Stopping on convergence

2. Application of Tabu search for optimal control of a polymerization reactor

The determination of open-loop time varying control policies that maximize or minimize a given performance index is referred as optimal control. The optimal control policies that guarantee the product property requirements and the operational constraints can be computed off-line, and are executed on-line in such a way that the process is operated in accordance with these control policies. Achievement of product quality is a major issue in polymerization processes since the molecular or morphological properties of a polymer product majorly influences its physical, chemical, thermal, rheological and mechanical properties as well as polymer applications. In a free radical copolymerization, composition drifts caused by variations in reactivities of comonomers can be mitigated by continuous addition of more reactive monomer while maintaining a constant mole ratio. The end use properties of the polymer product such as flexibility, strength and glass transition temperature are affected by the copolymer composition. Further, the molecular weight (MW) and molecular weight distribution (MWD) affects the important end use properties such as viscosity, elasticity, strength, toughness and solvent resistance. Hence, the determination of optimal control trajectories is important for the operation of polymerization reactors in order to produce a polymer with the desired product characteristics.

In the past, various methods have been reported for optimal control of polymerization reactors [12, 13, 14, 15, 16]. Most of the studies are based on classical methods of optimization such as Pontryagin’s maximum principle [17], which has been applied to solve the optimal control problems of different polymerization reactors [18, 19, 20, 21, 22, 23]. However, the classical methods have certain limitations for optimal control of polymerization reactors and these drawbacks have been discussed in literature [24]. The stochastic and evolutionary optimization methods are found beneficial over the conventional gradient-based search techniques because of their ability in locating the global optimum of multi-modal functions and searching design spaces with disjoint feasible regions.

2.1 Optimal control problem

The general open-loop optimal control problem of a lumped parameter batch/semi-batch process with fixed terminal time can be stated as follows. Find a control vector u(t) over t_f [t_o, t_f] to maximize (minimize) a performance index J(x,u):

$$J\left(x,u\right)=\underset{t_0}{\overset{t_f}{\int }}\phi \left[x\left(t\right),u\left(t\right),t\right]\mathit{dt}$$E4

Subject to

$$\dot{x}\left(t\right)=f\left(x\left(t\right),u\left(t\right),t\right),x\left(t_0\right)=x_0$$E5

$$h\left[x\left(t\right),u\left(t\right)\right]=0$$E6

$$g\left[x\left(t\right),u\left(t\right)\right]\le 0$$E7

$$x^L\le x\left(t\right)\le x^U$$E8

$$u^L\le u\left(t\right)\le u^U$$E9

In the above, J refers the performance index, x is the state variable vector, u is the control variable vector. Eq. (5) represents the system of ordinary differential equations with their initial conditions, Eqs. (6) and (7) specify the equality and inequality algebraic constraints and Eqs. (8) and (9) are the upper and lower bounds for the state and control variables.

2.2 Multistage dynamic optimization strategy

A multistage optimization results due to the natural extension of a single stage optimization. In a system involving multiple stages, the output from one stage becomes the input to the subsequent stage. The multistage dynamic optimization procedure can be referred elsewhere [24]. This type of optimization needs special techniques to split the problem into computationally manageable units. In multistage dynamic optimization, the optimal control problem of the entire batch duration is divided into finite number of time instants referred to as discrete stages. The control variables and the corresponding state variables that satisfy the objective function are evaluated in stage wise manner.

In this work, a multistage dynamic optimization strategy with sequential implementation procedure is presented for optimal control of polymerization reactors. The procedure for solving the optimal control problem is similar to that of the dynamic programming based on the principle of optimality [25]. The meta heuristic features of tabu search (TS) are exploited by implementing this strategy on a semi-batch styrene–acrylonitrile (SAN) copolymerization reactor. The procedure involves in discretizing the process into N stages, defining the objective function, f ⁱ, the control vector, uⁱ and the state vector, xⁱ for stage i. This procedure is briefed in the following steps:

1. The optimum value of the objective function, f¹[x¹] for stage 1 driven by the best control vector u¹ along with the state vector x¹ is represented as

$$\begin{array}{l}{f}^1\left[x^1\right]=min{f}_o^1\left[x^1,u^1\right]\\ u^1\end{array}$$E10

2. The value of the objective function, f²[x²] for stage 2 is determined based on the best control vector u² along with the state vector x² as given by

$$\begin{array}{l}{f}^2\left[x^2\right]=f^1\left[x^1\right]+min{f}_o^2\left[x^2,u^2\right]\\ u^2\end{array}$$E11

3. Recursive generalization of the above procedure for the k^th stage is represented by

$$\begin{array}{l}{f}^k\left[x^k\right]=f^{k-1}\left[x^{k-1}\right]+min{f}_o^k\left[x^k,u^k\right]\\ u^k\end{array}$$E12

In this procedure, $f_o^k$ represents the performance index of stage k.

2.3 The description of polymerization process and its mathematical representation

The solution copolymerization of styrene and acrylonitrile taking place in a semi-batch reactor is considered as a test bed for sequential implementation of the dynamic optimization strategy. The xylene is the solvent and AIBN is the initiator for the reaction. The feed which is a mixture of monomers, solvent, and initiator enters the reactor in semi-batch mode. The reactor initial volume is 1.01 L. The initially set design parameters are the solvent mole fraction f_s = 0.25 and initiator concentration I₀ = 0.05 mol/L. The mole ratio of monomers in the feed, M₁/M₂, is 1.5, where M₁ and M₂ are the molar concentrations of the unreacted monomers (styrene and acrylonitrile). The homogeneous solution free-radical copolymerization of styrene with acrylonitrile is described by the following reaction mechanism [26].

Initiation:

$$I\to 2R$$

$$R+M_1\stackrel{k_{i1}}{\to }{P}_{10}$$E13

$$R+M_2\stackrel{k_{i2}}{\to }{Q}_{01}$$

Propagation:

$$P_{n,m}+M_1\stackrel{k_{p11}}{\to }{P}_{n+1,m}$$

$$P_{n,m}+M_2\stackrel{k_{p12}}{\to }{Q}_{n,m+1}$$E14

$$Q_{n,m}+M_1\stackrel{k_{p21}}{\to }{P}_{n+1,m}$$

$$Q_{n,m}+M_2\stackrel{k_{p22}}{\to }{Q}_{n,m+1}$$

Combination termination:

$$P_{n,m}+P_{r,q}\stackrel{k_{\mathit{tc}11}}{\to }{M}_{n+r,m+q}$$

$$P_{n,m}+Q_{r,q}\stackrel{k_{\mathit{tc}12}}{\to }{M}_{n+r,m+q}$$E15

$$Q_{n,m}+Q_{r,q}\stackrel{k_{\mathit{tc}22}}{\to }{M}_{n+r,m+q}$$

Disproportionation termination:

$$P_{n,m}+P_{r,q}\stackrel{k_{\mathit{td}11}}{\to }{M}_{n,m}+M_{r,q}$$

$$P_{n,m}+Q_{r,q}\stackrel{k_{\mathit{td}12}}{\to }{M}_{n,m}+M_{r,q}$$E16

$$Q_{n,m}+Q_{r,q}\stackrel{k_{\mathit{td}22}}{\to }{M}_{n,m}+M_{r,q}$$

Chain Transfer:

$$P_{n,m}+M_1\stackrel{k_{f11}}{\to }{M}_{n,m}+P_{10}$$

$$P_{n,m}+M_2\stackrel{k_{f12}}{\to }{M}_{n,m}+Q_{01}$$

$$Q_{n,m}+M_1\stackrel{{k_f}_{21}}{\to }{M}_{n,m}+P_{10}$$E17

$$Q_{n,m}+M_2\stackrel{k_{f22}}{\to }{M}_{n,m}+Q_{01}$$

where P_n,m stand for a growing polymer chain with n units of monomer 1 (styrene) and m units of monomer 2 (acrylonitrile) with monomer 1 on the chain end. Similarly, Q_n,m refers to a growing copolymer chain with monomer 2 on the end. The M_n,m denotes inactive or dead polymer.

The molecular weight (MW) and molecular weight distribution (MWD) of the copolymer are computed using three leading moments of the total number average copolymers. The instantaneous k^th moment is given by:

$${{\lambda }_k}^d=\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }{\left({\mathit{nw}}_1+{\mathit{mw}}_2\right)}^k{M}_{n,m},k=0,1,2,.\dots \dots$$E18

where w₁ is the molecular weight of styrene and w₂ is the molecular weights of acrylonitrile. The total number average chain length (X_n), the total weight average chain length (X_w) and the polydispersity index (PD) are defined as:

$$X_n={{\lambda }_1}^d/{{\lambda }_0}^d$$

$$X_w={{\lambda }_2}^d/{{\lambda }_1}^d$$E19

$$\mathit{PD}=X_w/X_n$$

More details on modeling equations, reaction kinetics and numerical data pertaining to this system can be referred elsewhere [24].

The polymerization process model [24] is in the form of the general expression in Eq. (5). The set of state variables, X and the control vector U in the model are given by.

$$X\left(\mathrm{t}\right)={\left[M_1\left(t\right),M_2\left(t\right),I\left(t\right),V\left(t\right),{\lambda }_0\left(t\right),{\lambda }_1\left(t\right),{\lambda }_2\left(t\right)\right]}^{\mathrm{T}}$$

and

$$U\left(\mathrm{t}\right)={\left[T\left(t\right),u\left(t\right)\right]}^{\mathrm{T}}$$E20

In the above, I is the concentration of the unreacted initiator, V is the reaction mixture volume, λ_k (k = 0, 1, …) is the k^th moment of the dead copolymer molecular weight distribution, T is the reaction mixture temperature, an u is the volumetric flow rate of the feed mixture to the reactor.

2.4 Control objectives

The desired values of copolymer composition (F₁D) and number average molecular weight (MWD) are chosen to be 0.58 and 30000, respectively. Minimizing the deviations of copolymer composition, F₁ and molecular weight, MW from their respective desired values during the entire span of the reaction are specified as the desired objectives. In order to attain these objectives, the control variables are set as monomer addition rate, u and reactor temperature, T. If one manipulative variable is used to control one polymer quality parameter, the uncontrolled property parameters may deviate from their desired values as the reaction proceeds. The optimal control problem involves in optimizing the single objectives as well as both the objectives simultaneously. The objectives of SAN copolymerization are specified as.

$$J_1={\left[1-\mathit{MW}\left(t\right)/\mathsf{MWD}\right]}^2$$E21

$$J_2={\left[1-F_1\left(t\right)/F_{1}D\right]}^2$$E22

$$J_3={\left[1-F_1\left(t\right)/F_{1}D\right]}^2+{\left[1-\mathit{MW}\left(\mathrm{t}\right)/\mathsf{MWD}\right]}^2$$E23

Here MW and F₁ are the molecular weight and copolymer composition, and MWD and F₁D are their respective desired values. The notation t here refers discrete time. The hard constraints are set as.

$$0\le u\left(t\right)\le 0.07\left(\mathrm{l}/min\right)$$

$$320\le T\left(t\right)\le 368\left(k\right)$$E24

$$V\left(t\right)\le 4.0\mathrm{l}$$

These constraints on operating variables are selected by taking into consideration of reaction rate, heat transfer limitation and reactor safety. Determination of temperature policy, T, to satisfy J₁ (Eq. (21)) and monomer feed policy, u, to satisfy J₂ (Eq. (22)) are considered as single objective optimization problems. Determination of both u and T policies that satisfy J₃ (Eq. (23)) is considered as a problem of simultaneous optimization.

2.5 Design and implementation

The design and implementation of tabu search for optimal control of copolymerization reactor is explained as follows. The total span of reaction time in SAN copolymerization reactor is split into finite number of time instants referred to discrete stages. The total time of reaction is fixed at 300 min. The duration of reaction time is divided into 19 stages with 20 time points, each stage having a duration of 15 min. Thus the discrete control sequences of feed flow rate, u and temperature, T are specified as.

$$u={\left[u_1,u_2,..\dots ,u_{20}\right]}^{\mathrm{T}}$$E25

$$T={\left[T_1,T_2,..\dots ,T_{20}\right]}^{\mathrm{T}}$$E26

The control sequences are located at equal distances. Initially, the control vector is specified as constant value at each of the discrete stages. The control input at the beginning of the first stage is chosen as the lower bound of the input space. The elements of tabu search for computing the optimal control policies are specified as follows. The sizes of both recency and frequency based tabu lists are set as 50. The sizes of these lists are chosen such that they prohibit revisiting of un-prosperous solutions in the search process. An intensification procedure in the form of a sine function is used. The oscillation period α of the sine function is controlled by a parameter (θ) whose value is 4.0001. A sigmoid function based aspiration criterion with the parameters as k_center = 0.3 and σ = 7/M with M as specified number of iterations is used. For optimizing MW, ten neighbors are formed at the end of the first stage by considering random changes in the search space of T with an incremental change of −0.4 to 0.4. The number of neighbors specified for each control input of each stage is 10. The integration of process model is performed for a time period of 15 min with a time step of 1 min from the starting to the end of the first stage and the objective function values are calculated for all the generated neighbors at the end of this stage. The iterative convergence of TS leads to establish the best control input (T) along with its objective function. This provides the optimal control point for the first stage, which is then used as a starting point for the second stage solution. For second stage solution, random neighbors are created at the end of the second stage around the optimal T of the first stage. The integration of the model is performed from the beginning to the end of the first stage based on the initial control point (T) and from the starting to the end of the second stage for each of the neighbors generated at the end of second stage. The optimal control inputs for successive stages are computed until the end of last stage in a similar manner. The control input values that are computed at the end of each stage thus represents the optimal control policy for T. For F₁ optimization, ten neighbors are generated at the end of first stage with an incremental variation of −1.0 x 10⁻⁶ to 1.0 x 10⁻⁶. In analogous manner, optimal control policy for u(F₁) is found by following the similar TS procedure as in T policy. For determining the dual control policies of T and u, multistage dynamic optimization by TS is performed by accounting incremental variations in neighbors generation of T and u within the ranges of −0.4 to 0.4, and − 1.0 x 10⁻⁶ to 1.0 x 10⁻⁶, respectively. This case of multistage dynamic optimization involves the computation of the objective function values for 100 neighbor combinations at each of the control point corresponding to T and u.Figure 2 shows the implementation of TS strategy for optimal control of SAN copolymerization reactor.

Figure 2.

Multistage dynamic optimization of copolymerization reactor using tabu search.

2.6 Analysis of results

Tabu search (TS) is designed and applied to compute the optimal control policies that meet the single and multiple objectives of SAN copolymerization reactor. The dual control policies of T and u computed by TS for multistage dynamic optimization of polymerization reactor along with the objective function values are shown in Figure 3. These results exhibit the efficiency of TS in computing both the temperature and feed rate policies for the reactor to maintain MW and F₁ almost near to their desired values. The computational efficiencies of TS are evaluated by means of execution times, normalized absolute error values and memory storage requirements as shown in Table 1.

Figure 3.

Dual control policies and objectives: (a) Temperature policy, (b) Feed rate policy, (c) Molecular weight as an objective, and (d) Copolymer composition as an objective.

Strategy	Control Policy	Computational efficiency		Normalized absolute error	Memory storage, k
		Convergence time (sec)	Implementation time (sec)
TS	T	1.89	16.83	0.000566	2224
	U	2.19	13.81	0.01379	2076
	T & u	4.25	19.21	0.007535	3084

Table 1.

Computational efficiencies of TS strategy.

3. Application of Tabu search for inverse modeling of anaerobic fixed bed biofilm reactor

Tabu search is applied to determine the parameters of kinetic and film thickness models through the validation of the mathematical models of the process with the help of measured data acquired from an experimental fixed bed anaerobic biofilm reactor used in the treatment of pharmaceutical industry wastewater.

3.1 Experimental biofilm reactor and its description

The schematic of a laboratory scale fixed bed biofilm reactor setup used for industry waste water treatment is shown in Figure 4. The description of the experimental unit, its auxiliary items and the packing specifications are given as follows.

Figure 4.

Experimental biofilm reactor setup.

The reactor consists of a QVF glass column of 1 m height and 0.1016 m internal diameter. A Teflon perforated plate of thickness 3 mm is provided at the base of the column to support the packing material. Two openings fitted with valves made of Teflon are provided at the bottom of the reactor so that one valve is used to control the flow rate of influent pumped into the reactor and the other valve is used to discharge the excess sludge accumulated. The top of the column has a provision to place a thermometer to measure the temperature. Three sampling ports are placed at equal distances along the side of the column to withdraw samples. The temperature of the column is maintained by varying the voltage through the 0.2 kW heating tape connected to a dimmer-stat. The column is insulated with 1.0 in. asbestos rope. The entire setup is supported by a 1.0 in. M.S. pipeline structure. The column and the packing specifications are given elsewhere [27].

3.2 Experiments and data generation

The experimental setup shown in Figure 4 is used to conduct experiments for anaerobic treatment of pharmaceutical industry wastewater involving different packing materials with varying organic loading rates and feed rates, and different hydraulic retention times. The influent kept in a 2 l capacity vessel is fed to the column through an inlet connection located at the bottom using a peristaltic pump made of Watson Marlow. The effluent sample is collected from an outlet located at the top of the column. The column and the packing specifications are given elsewhere [27]. The reactor is filled with the packing material and then seeded with 1.0 l of active anaerobic sludge brought from M/s. Alkabir, Hyderabad, India. To maintain the initial growth of microorganisms, the reactor is further added with 2.6 l of very dilute solution of synthetic glucose medium (1000 mg/l) with nutrients such as total nitrogen as urea (125 mg/l), total phosphorous as KH2PO4 (50 mg/l), NaCl (50 mg/l), KCl (40 mg/l), CaCl2 (15 mg/l), MgCl2 (10 mg/l) and FeCl3 (2 mg/l). The pH of the solution is maintained at 7.2 by adding 0.1 N NaOH. The reactor is made biologically active by keeping it undisturbed for about 20 days. The gas (CH4 + CO2) evolved at the top of the column confirms the biological activity of microorganisms.

The wastewater procured from a typical pharmaceutical industry is a complex medium and its composition is reported elsewhere [28]. Substrates of varying COD concentrations are prepared using the wastewater in order to use them as feed solutions for biofilm reactor experiments. Once the biological activity of the bed is detected, the diluted industry wastewater solution equivalent to three bed volumes is sent to the reactor by an electronic dosing pump. NaHCO₃ is added to maintain the pH of the feeding solution at 7.2. Different experiments are conducted for treating the pharmaceutical industry wastewater using the feed solutions having the COD concentrations of 4,980 mg/l, 11,860 mg/l, 23,460 mg/l and 40,720 mg/l. The hydraulic retention time (HRT) of each experiment is varied from 2 to 12 days based on the substrate concentration of the feed solution. Each experimental run subsequent to attainment of its HRT is allowed to continue for 13 days so as the operation reaches stable state condition by that time.

Thirteen experiments are conducted using the different feed solutions prepared from the pharmaceutical industry wastewater. The reactor is provided with the sampling ports to withdraw samples at bed heights of 0.25 m, 0.5 m and 0.75 m from the bottom. Samples collected from the three sampling ports of the column as well as from the effluent water are analyzed for COD, BOD, TVA, TA, pH etc. The stable state condition of each experiment is observed with respect to the minimal variation in reactor effluent pH, alkalinity and COD concentrations. Thus 13 experiments are conducted under different feed stream conditions. The data corresponding to the influent and effluent COD concentrations, HRT and OLR of all experiments are given elsewhere [28].

3.3 Mathematical and kinetic models

The application of mathematical models to biofilm reactors suffers due to lack of availability accurate kinetic models and uncertainty in model parameters. Therefore, appropriate selection of kinetic model and accurate determination of its parameters is required for successful modeling of biofilm reactors. The biological reaction rates are generally represented by Monod kinetics and also occasionally by Contois and Haldane models. In the past, efforts have been made to determine the parameters of kinetic models of biofilm reactors, mostly experimentally [29, 30, 31]. Various researchers [32, 33, 34, 35] have reported that numerical evaluation of kinetic parameters is an attractive alternative to the experimental methods for biofilm processes. By numerical approach, the parameters of the kinetic models of biofilm processes are determined as a consequence of the validation of their mathematical models with the aid of measured data. The approach of determining the model parameters such that the behavior of the process model approximates the observed process behavior is known as inverse modeling.

3.3.1 Mathematical models

Mathematical models of varying complexities involving different types of kinetic and film thickness models are used to represent the wastewater treatment process in the biofilm reactor.

3.3.1.1 One dimensional model

This model accounts the rate of mass transfer to be proportional to the concentration difference between the interface and the bulk fluid. This model assumes no accumulation of the component at the interface under steady state condition. The differential equation governing the fluid phase is given by.

$$-u\frac{\mathit{dc}}{\mathit{dz}}=k_g{a}_v\left(c-c_s^s\right)$$E27

with the boundary condition:

$$\mathrm{c}={\mathrm{c}}_0\mathrm{at}\mathrm{z}=0$$E28

In this model, the substrate concentration in the bulk fluid varies only in the axial direction.

The differential equation governing the solid phase is given by.

$$\frac{D_f}{{\xi }^{a-1}}\frac{d}{\mathit{d\xi }}\left({\xi }^{a-1}\frac{{\mathit{dc}}_s}{\mathit{d\xi }}\right)-r_s\left(c_s\right)=0$$E29

with the boundary conditions:

$$\frac{{\mathit{dc}}_s}{\mathit{d\xi }}=0\mathrm{at}\xi =0$$

$$k_g\left(c_s^s-c\right)=-D_f\frac{{\mathrm{dc}}_{\mathrm{s}}}{\mathrm{d}\xi }\mathrm{at}\xi =L_f$$E30

The superscript ‘a’ in Eq. (29) refers to 1, 2, and 3 for planar, cylindrical and spherical geometries, respectively, and D_f denotes substrate diffusion coefficient in the biofilm (m²/day). The notation for other terms in the above equations is denoted as k_g is mass transfer coefficient (m/s), a_v is specific surface area of particle (m⁻¹), c is substrate concentration in the bulk fluid (kg/m³), c_s is substrate concentration in the bio-film (kg/m³), c^s_s is substrate concentration on the bio-film surface (kg/m³), ξ is space coordinate in the bio-film (m), z is axial coordinate (m) and L_f is thickness of the bio-film (m).

3.3.1.2 Two dimensional model

This model considers the variation of the substrate concentration in bulk fluid in both the axial and radial directions. The substrate transfer occurring through the biofilm is characterized by physical diffusion and reaction, and this process is described by second order differential equation. Pressure losses along the length of the reactor are neglected. The differential equation governing the bulk fluid phase is given by.

$$u\frac{\mathit{dc}}{\mathit{dz}}=\frac{\mathit{\epsilon D}}{r}\left[\frac{\partial }{\partial r}\left(r\frac{\partial c}{\partial r}\right)\right]-k_g{a}_v\left(c-c_s^s\right)$$E31

where D is substrate diffusion coefficient in the bulk fluid (m²/day), u is superficial velocity (m/s), r is radial coordinate (m) and ε is porosity of bed. The boundary condition for solving this equation is given in Eq. (28). The solid phase biofilm model for two dimensional model is same as in one dimensional model as in Eq. (29) and its boundary conditions are expressed by Eq. (30).

3.3.1.3 Mass transfer coefficient

The substrate mass transport occurring from the bulk fluid to the biofilm surface across the diffusion layer is defined by

$$j_D=\frac{k_g{\left(\mathit{sc}\right)}^{2/3}}{u}$$E32

where Sc is Schmidt number (μ/ρD). The j_D-factor is calculated using the following correlation:

$$j_D\epsilon =\frac{0.765}{{Re}^{0.82}}+\frac{0.365}{{Re}^{0.386}}$$E33

$$\text{where}Re=\frac{\rho {\mathit{ud}}_{\mathrm{p}}}{\mu }\text{for}0.01<Re<15,000$$

where d_p is equivalent particle diameter (m). The mass transfer coefficient, k_g obtained from Eq. (32) is used in bulk fluid phase equations, Eq. (27) and Eq. (31).

3.3.2 Kinetic models

Various kinetic models can be used to represent the bioprocess kinetics [28], of which the Haldane and Edward models are given as follows.

3.3.2.1 Haldane model

This rate expression is generally valid when the substrate concentration is high. According to this model, the substrate consumption rate, r_s in biofilm is given by the following equation:

$$r_s=\left({\mu }_{\max }{\rho }_s/Y\right)\frac{C_s}{K_s+C_s+\frac{C_s^2}{K_I}}$$E34

where C_s is substrate concentration, μ_max is maximum specific growth rate, ρ_s is density of biomass, Y is yield coefficient, K_s is half velocity constant, and K_I substrate inhibition constant.

3.3.2.2 Edward model

This model considers substrate inhibition, according to which the substrate consumption rate, r_s in biofilm is given by the equation:

$$r_s=\left(\frac{{\mu }_{\max }{\rho }_s}{\mathrm{Y}}\right)\left(exp\left(\frac{-C_s}{K_I}\right)-exp\left(\frac{-C_s}{K_s}\right)\right)$$E35

3.3.2.3 Film thickness

The biofilm thickness (L_f) is has significant influence on substrate conversion. A uniform film thickness with fixed value cannot represent the realistic situation and the film thickness L_f is expected to vary along the reactor [35]. The varying thickness of the biofilm is determined by the substrate flux from the bulk fluid into the biofilm as wellas the growth and decay rates of the bacteria. It is assumed that the biofilm is composed of a static portion and a variable portion. The variable portion is supposed to raise with the organic loading rate and lessen with the hydraulic loading rate. The above point of view leads to the following relations for estimating the biofilm thickness:

$$L_f=a+b\mathit{OLR}$$E36

$$L_f=a+b\mathit{OLR}-\frac{c}{\mathit{HLR}}$$E37

where OLR is organic loading rate (kg/m³/day), HLR hydraulic loading rate (m/day), and a, b and c are constants to be determined.

3.4 Procedure for solving the model equations

To facilitate the solution of mathematical models, the height of the column (1 m) is divided into 50 equal steps, each step representing 0.02 m. In one dimensional model, the fluid phase equation (Eq. (27)) along with its boundary condition (Eq. (28)) is solved for each height step in the axial direction by using Runge–Kutta (RK) 4th order method. The solid phase equation for the biofilm (Eq. (29)) is solved for each segment using orthogonal collocation on finite elements (OCFE) [36]. Two finite elements, each with four internal collocation points are considered for OCFE implementation. In two dimensional model, the fluid phase equation (Eq. (31)) is transformed to ODE by finite difference technique and the resulting equation along with its boundary condition is solved by using 4th order RK method for each height step of the axial direction. The solution for solid phase equation of this model is same as in one dimensional model.

3.5 Tabu search for inverse modeling of biofilm reactor

The implementation procedure of tabu search is given in Section 1.3. In this work, tabu search is used to determine the parameters of kinetic and film thickness models in association with the validation of the mathematical models using the data of measurements obtained from an experimental fixed bed anaerobic biofilm reactor. The parameter estimation via inverse modeling is carried out by defining an objective function, J given by.

$$J=f\left(\alpha \right)=\sum _{i=1}^l{\left(y_i-{\widehat{y}}_i\right)}^2$$E38

where α is the vector of parameters, l is the number of observations, y_i is the measured value of the i^th variable, and ${\widehat{y}}_i$ is the corresponding predicted value. Tabu search is applied to estimate the parameters of kinetic and film thickness expressions of different modeling configurations with the support of one dimensional (1D) and two dimensional (2D) mathematical models of the biofilm reactor. Random variation in all the parameters is considered to generate neighbors so as to provide the maximum improvement to the current solution. Recency based tabu list with a length of 100 and frequency based tabu list with a length of 50 are employed. A sigmoid function based aspiration criterion, Eq. (2) is employed with k_center as 0.3 and $\sigma$ as 7/M, where M refers specified number of iterations. An intensification strategy with the coefficient α taking the format of a sine function, Eq. (1) is employed, in which the value of θ is assigned to be 4.0001. The maximum number of iterations are considered as 100. The termination-on-convergence criteria, Eq. (3) is used as stopping criterion.

3.6 Analysis of results

Parameter estimation by tabu search is performed by devising different modeling configurations to represent the biofilm reactor. These configurations involve the Edward kinetics and Haldane kinetics with substrate inhibition along with the two and three parameter film thickness expressions. The modeling configuration that is formed by 1D model along with Haldane kinetics and two parameter film thickness expression is referred to 1D-Haldane-2film. Other modeling configurations are formed and abbreviated in the same manner. The effectiveness of these modeling configurations are assessed by comparing the model predictions with experimental data and further by means of performance measures such as cost function (CF) and model efficiency (ME) [28]. A conventional optimization method called Nelder–Mead optimization (NMO) [37, 38] is also employed for comparison with the tabu search. The results are analyzed to find the better suitability of mathematical models using the quantitative performance measures (CF and ME). These results evaluated for 1D and 2D models with Haldane and Edward kinetics involving two film and three film expressions indicate the better suitability of 2D model over 1D model. When the results are analyzed to assess the better suitability optimization algorithm, the analysis of results of different modeling configurations indicate the effectiveness of TS over NMO. The results evaluated to find the usefulness of kinetic models indicate the better suitability of Edward kinetics over Haldane kinetics. The results analyzed to assess the film thickness expressions indicate the appropriateness of the three parameter film thickness expression over two parameter expression. The iterative convergence of parameters estimated by TS and NMO are shown in Figure 5. The comparison of the prediction results of 2D model with Edward kinetics and three parameter film thickness expression evaluated using TS and NMO with those of experimental results are shown in Figure 6. The analysis of the results show the better performance of TS over NMO for inverse modeling of biofilm reactor. From these results, it is found that the TS involving 2D model, Edward kinetics and three parameter film thickness expression better represents the fixed bed biofilm reactor involved in the treatment of pharmaceutical industry wastewater.

Figure 5.

Iterative convergence of parameter estimates by TS and NMO: (a) kinetic parameters, (b) film thickness parameters.

Figure 6.

Comparison of experimental results with model predicted substrate conversions of 2D model with Edward kinetics.

4. Conclusions

Tabu search is a stochastic optimization method used to solve global optimization problems. The effectiveness of tabu search for modeling and optimization is explored with the applications concerning to chemical and environmental systems of varying complexities. The meta heuristic Tabu search (TS) is designed and applied to determine the optimal control policies of a SAN copolymerization reactor and inverse modeling of a biofilm reactor. The computational efficiencies of TS are evaluated in terms of execution times, normalized absolute error values and memory storage requirements. The results of polymerization reactor show the usefulness of TS in determining the optimal temperature and feed rate policies to maintain the objectives MW and F₁ near to their desired values. The results of biofilm reactor explain the efficacy of TS in appropriately configuring the kinetic and film thickness models as a consequence of validation of mathematical models.