Parameter Identification of Anaerobic Wastewater Treatment Bioprocesses Using Particle Swarm Optimization

Abstract

This paper deals with the offline parameters identification for a class of wastewater treatment bioprocesses using particle swarm optimization (PSO) techniques. Particle swarm optimization is a relatively new heuristic method that has produced promising results for solving complex optimization problems. In this paper one uses some variants of the PSO algorithm for parameter estimation of an anaerobic wastewater treatment process that is a complex biotechnological system. The identification scheme is based on a multimodal numerical optimization problem with high dimension. The performances of the method are analyzed by numerical simulations.

1. Introduction

It is well known that the biotechnology is one of the fields that over the last decades have a high development. Therefore, due to its advantages, the control of industrial bioprocesses has been an important practical problem attracting wide attention. Biotechnology applications can be found especially in agriculture, in food industry, in medicine and pharmaceutical processes, in waste treatment processes, and so forth. A frequent and important challenge in control of such living processes is finding an accurate model of the system. The bioprocesses are highly nonlinear, and their kinetic parameters are usually badly or inadequately known [1]. This problem becomes of great importance in complex systems where critical instability of the process must be avoided. Parameters characterizing the internal behavior of biotechnological systems are usually not directly accessible to measurement. Their measurement is usually approached indirectly as a parameter estimation problem [2]. In this paper a dynamic model describing the internal structure of the system is formulated, and an algorithm based on PSO for parameter estimation is designed.

In recent years, a progress has been made in the area of continuous-time system identification [3]. Even if the most physical systems are naturally continuous, a much more attention has been paid to parameter estimation of discrete-time systems, mainly because they are better suited for numerical implementations. Continuous-time identification makes possible a more direct link to the physical properties and operation of the underlying systems and the direct estimation of physical parameters which have a clear significance. The most common approach for parameter estimation of linear or nonlinear systems is the use of prediction-error identification methods (PEM) [4]. In this category falls the well-known least squares methods or the maximum likelihood methods. In this approach, identification consists in minimization of a scalar-valued function of the model parameter. In general, this function cannot be minimized by analytical methods so the solution has to be found by iterative, numerical techniques. There is an extensive literature on such numerical problems. In classical approach the most used procedures are the quasi-Newton methods and interior point algorithms. The main drawback of these nonlinear parameter optimization techniques is that they are often unreliable; for example, they give no guarantee of converging to a true minimum. The increasing computational power of personal computers and microcontrollers allowed the implementation of several optimization algorithms inspired from natural phenomena. Examples of these algorithms include the Simulated Annealing [5], Genetics Algorithms (GA) [6], or Ant Colony Optimization [7] algorithms. Particle Swarm Optimization (PSO) [8] is among these nature inspired algorithms. It is inspired by the ability of birds flocking to find food that they have no previous knowledge of its location. Every member of the swarm is affected by its own experience and its neighbors’ experiences. Although the idea behind PSO is simple and can be implemented by two lines of programming code, the emergent behavior is complex and hard to completely understand [9, 10].

The most important approaches for the yield and kinetic coefficients estimation of biotechnological systems make use of the state transformations based on the general structure [11]. In this paper we propose an identification method based on particle swarm optimization techniques for these classes of biotechnological systems considering that the unknown parameters can appear in rational relations with measured variables. The paper is organized in the following way. The nonlinear dynamical model of an anaerobic wastewater treatment bioprocess is given in Section 2. Section 3 presents the identification algorithm using the particle swarm optimization techniques. Some numerical simulations are presented in Section 4 and conclusions in Section 5.

2. Nonlinear Dynamical Model of Anaerobic Wastewater Treatment Bioprocesses

A process that takes place in a bioreactor can be described as a set of 𝑚 biochemical reactions involving 𝑛 components (with 𝑛≥𝑚). The global dynamics can be represented by the following dynamical state-space model [1]: 𝑑𝜉()𝑑𝑡=𝐾⋅𝜑𝜉,𝑡−𝐷𝜉+𝐹−𝑄, (1) where 𝜉∈ℜ^𝑛×1.

This model describes the behavior of an entire class of biotechnological processes and is referred to as the general dynamical state-space model of this class of bioprocesses [9]. In (1), the term 𝐾⋅𝜑(𝜉,𝑡) is the rate of consumption and/or production of the components in the reactor; that is, the reaction kinetics and the term −𝐷𝜉+𝐹−𝑄 represents the exchange with the environment. The strongly nonlinear character of the model (1) is given by the reaction kinetics. In many situations, the yield coefficients, the structure, and the parameters of the reaction rates are partially known or unknown. Many of the evolved control methods for these kinds of systems—like model predictive control and robust or adaptive control—are based on good initial estimates of the yield and kinetic parameters.

Anaerobic digestion is a multistage biological wastewater treatment process whereby bacteria, in the absence of oxygen, decompose organic matter to carbon dioxide CO₂, methane CH₄, and water [12]. Four metabolic paths can be identified in this process: two for acidogenesis and two for methanation (see Figure 1). In the first acidogenic path, glucose (or another complex substrate) from the wastewater is decomposed into volatile fatty acids (acetates, propionic acid), hydrogen H₂, and inorganic carbon by acidogenic bacteria.

Figure 1 

A schematic view of a typical anaerobic digestion process.

In the second acidogenic path, Obligate Hydrogen Producing Acetogens (OHPA) decompose propionate into acetate, H₂, and inorganic carbon. In the first methanation path, acetate is transformed into CH₄ and CO₂ by acetoclastic methanogenic bacteria, while, in the second methanation path, H₂ combines inorganic carbon to produce CH₄ under the action of hydrogenophilic methanogenic bacteria. The reaction scheme of this complex bioprocess involves 4 reactions and 10 components.

Since the anaerobic digestion is a complex bioprocess, this dynamical model being described by ten differential equations, for control purpose an appropriately reduced order model can be used. Using the singular perturbation method, the following reduced order dynamical model can be obtained: 𝑋̇₁=𝜑₁−𝐷𝑋₁,𝑆̇₁=−𝑘₁𝜑₁−𝐷𝑆₁+𝐷𝑆_in,𝑋̇₂=𝜑₂−𝐷𝑋₂,𝑆̇₂=𝑘₃𝜑₁−𝑘₂𝜑₂−𝐷𝑆₂,𝑆̇₅=𝑘₄𝜑₁+𝑘₅𝜑₂−𝐷𝑆₅−𝑄_CO2,𝑃̇=𝑘₆𝜑₂+𝑘₇𝜑₁−𝐷𝑃−𝑄_𝑃, (2)(3)(4)(5)(6)(7) where 𝑋₁, 𝑋₂ are acidogenic and acetoclastic methanogenic bacteria, respectively, and 𝑆₁, 𝑆₂, 𝑆₅ are glucose, acetate and inorganic carbon, respectively, and 𝑃 is methane; 𝜑₁, 𝜑₂ are the rates of first acidogenic reaction and methanation reaction respectively, 𝑄_𝑃=𝑐_𝑝𝑃 with 𝑐_𝑝>0 and 𝑄_CO2 represent gaseous outflow rates of CH₄ and CO₂, respectively, 𝑆_in is the influent substrate concentration, 𝐷 the dilution rate, and 𝑘_𝑖(𝑖=0,…,7) are yield coefficients. Each reaction rate is a growth rate and may be written as 𝜑_𝑖=𝜇_𝑖𝑋_𝑖,𝑖=1,2, where 𝜇_𝑖,𝑖=1,2, is the specific growth rate of reaction 𝑖. Defining the state vector as 𝜉=𝑋₁𝑆₁𝑋₂𝑆₂𝑆₅𝑃^𝑇, the model (2)–(7) can be written in matrix form as 𝑋̇₁𝑆̇₁𝑋̇₂𝑆̇₂𝑆̇₅𝑃̇=10−𝑘₁0𝑘01₂−𝑘₃𝑘₄𝑘₆𝑘₅𝑘₇𝜑₁𝜑₂𝑋−𝐷¹𝑆₁𝑋₂𝑆₂𝑆₅𝑃+0𝐷𝑆_in0000−0000𝑄_CO2𝑄_𝑃, (8) or in compact form as 𝜉̇=𝐾𝜑(𝜉)−𝐷𝜉+𝐹−𝑄, (9) where 𝐹=0𝐷𝑆_in0000^𝑇 is the vector of inflow rates, 𝑄=0000𝑄_CO2𝑄_𝑃^𝑇 is the vector of gaseous outflow rates, 𝜑=𝜑₁𝜑₂^𝑇 is the vector of reaction rates, which can be written as 𝜑(𝜉)=𝐺(𝜉)𝛼(𝜉), with 𝐺(𝜉) being a diagonal matrix whose entries are products of the component concentrations involved in each reaction and 𝛼=𝛼₁𝛼₂^𝑇 the vector of specific reaction rates, and 𝐾 is the yield coefficient’s matrix. The matrices 𝐾 and 𝐺 have the following structure: 𝐾=1−𝑘₁0𝑘₂𝑘₄𝑘₆001−𝑘₃𝑘₅𝑘₇^𝑇,𝑋𝐺=¹𝑆₁00𝑋₂𝑆₂. (10)

The most difficult task for the construction of the dynamical model is the modeling of the reaction kinetics [13]. The form of kinetics is complex, nonlinear, and in many cases unknown. In our study one considers that reaction rates are given by the Monod law 𝜑₁(𝜉)=𝜇^∗₁𝑆₁⋅𝑋₁𝐾_𝑀1+𝑆₁, (11) and the Haldane kinetic model 𝜑₂(𝜉)=𝜇^∗₂𝑆₂⋅𝑋₂𝐾_𝑀2+𝑆₂+𝑆²₂/𝐾_𝑖, (12) where 𝐾_𝑀1,𝐾_𝑀2 are Michaelis-Menten constants, 𝜇^∗₁,𝜇^∗₂ represent specific growth rate coefficients, and 𝐾_𝑖 is the inhibition constant.

For simplicity, we will denote the unknown plant parameters by the vector 𝜃=𝜃₁𝜃₂⋯𝜃₁₂^𝑇, (13) where 𝜃₁=𝑘₁;𝜃₂=𝑘₂;𝜃₃=𝑘₃;𝜃₄=𝑘₄;𝜃₅=𝑘₅𝜃₆=𝑘₆𝜃₇=𝑘₇𝜃₈=𝜇^∗₁;𝜃₉=𝜇^∗₂;𝜃₁₀=𝐾_𝑀1;𝜃₁₁=𝐾_𝑀2;𝜃₁₂=𝐾_𝑖. (14)

3. Parameter Estimation Using PSO

At the beginning of parameter estimation, the input and output data are known and the real system parameters are assumed as unknown. The identification problem is formulated in terms of an optimization problem in which the error between an actual physical measured response of the system and the simulated response of a parameterized model is minimized. The estimation of the system parameters is achieved as a result of minimizing the error function by the PSO algorithm.

3.1. Problem Statement

Consider the following n-dimensional nonlinear system: 𝑑𝜉(𝑡)(),𝑑𝑡=𝑓𝜉,𝑡;𝜃 (15) where 𝜉∈𝑅^𝑛 is the state vector, 𝜃∈𝑅^𝑚 is the unknown parameters vector, and 𝑓 is a given nonlinear vector function.

To estimate the unknown parameters in (15), a parameter identification system is defined as follows: 𝑑𝜉̂(𝑡)̂𝜃,𝑑𝑡=𝑓𝜉,𝑡; (16) where 𝜉̂∈𝑅^𝑛 is the estimated state vector and 𝜃∈𝑅^𝑚 is the estimated parameters vector.

The objective function defined as the mean squared errors between real and estimated responses for a number 𝑁 of given samples is considered as fitness of estimated model parameters: 1𝑉=𝑁+𝐷^𝐷𝑁_𝑗=1𝑘=1𝜉^𝑘_𝑗−𝜉̂^𝑘_𝑗², (17) where 𝐷 is the number of measurable states, 𝑁 is the data length used for parameter identification, whereas 𝜉^𝑘_𝑗 and 𝜉̂^𝑘_𝑗 are the real and estimated values of state 𝑗 at time 𝑘, respectively.

This objective function is a difficult function to minimize because there are many local minima, and the global minimum has a very narrow domain of attraction. Our goal is to determine the system parameters, using particle swarm optimization algorithms in such a way that the value of 𝑉 is minimized, approaching zero as much as possible.

3.2. Overview of Basic PSO Algorithms

During the last decade, PSO algorithms have gained much attention and wide applications in different fields due to their effectiveness in performing difficult optimization issues, as well as simplicity of implementation and ability to fast converge to a reasonably good solution [14–16]. PSO is a population-based heuristic global optimization technique, first introduced by Kennedy and Eberhart [8] and referred to as a swarm-intelligence technique. It is motivated from the simulation of social behavior of animals such as bird flocking, fish schooling, and swarm. In this algorithm, the population is called a swarm, and the trajectory of each particle in the search space is controlled through the medium of a term called “velocity,” according to its own flying experience and swarm experience in the search space. Mathematical description of basic PSO and some important variants is presented in the following.

Candidate solutions of a population called particles coexist and evolve simultaneously based on knowledge sharing with neighboring particles. Each particle represents a potential solution to the optimization problem, and it has a fitness value decided by optimal function. Supposing that search space is 𝐷-dimensional, each individual is treated as a particle in the 𝐷-dimensional search space. The position and rate of position change for 𝑖th particle can be represented by 𝐷-dimensional vector, 𝑥_𝑖=(𝑥_𝑖1,𝑥_𝑖2,…,𝑥_𝑖𝐷) and 𝑣_𝑖=(𝑣_𝑖1,𝑣_𝑖2,…,𝑣_𝑖𝐷), respectively. The best position previously visited by the 𝑖th particle is recorded and represented as 𝑝_𝑖=(𝑝_𝑖1,𝑝_𝑖2,…,𝑝_𝑖𝐷), called pbest. The swarm best position previously visited by all the particles in the population is represented as 𝑝_𝑔=(𝑝_𝑔1,𝑝_𝑔2,…,𝑝_𝑔𝐷), called gbest. Then particles search their best position, which are guided by swarm information 𝑝_𝑔 and their own information 𝑝_𝑖. Each particle modifies its velocity to find a better solution (position) by applying its own flying experience (i.e., memory of the best position found in earlier flights) and the experience of neighboring particles (i.e., the best solution found by the population). Each particle position is evaluated by using fitness function and updates its position and velocity according to the following equations: 𝑣_𝑖^𝑘+1=𝜔⋅𝑣^𝑘_𝑖+𝑐₁𝑟₁𝑝𝑏𝑒𝑠𝑡^𝑘_𝑖−𝑥^𝑡_𝑖+𝑐₂𝑟₂𝑔𝑏𝑒𝑠𝑡^𝑘_𝑖−𝑥^𝑘_𝑖,𝑥_𝑖^𝑘+1=𝑥^𝑘_𝑖+𝑣_𝑖^𝑘+1, (18) where 𝑘 is iteration number, 𝜔 is inertia weight, 𝑐₁ and 𝑐₂ are two acceleration coefficients regulating the relative velocity toward local and global best positions, and 𝑟₁ and 𝑟₂ are two random numbers from interval [0,1]. Many effects have been made over the last decade to determinate the inertia weight. Various studies has shown that under certain conditions convergence is guaranteed to a stable equilibrium point. These conditions include 𝜔>(𝑐₁+𝑐₂)/2−1 and 0<𝜔<1. The technique originally proposed was to bound velocities so that each component of 𝑣_𝑖 is kept within the range [−𝑉_max,+𝑉_max].

Unfortunately, this simple form of PSO suffers from the premature convergence problem, which is particularly true in complex problems since the interacted information among particles in PSO is too simple to encourage a global search. Many efforts have been made to avoid the premature convergence. One solution is the use of a constriction factor to insure convergence of the PSO, introduced in [17]. Thus, the expression for velocity has been modified as 𝑣_𝑖^𝑘+1=ℎ⋅𝑣^𝑘_𝑖+𝑐₁𝑟₁𝑝𝑏𝑒𝑠𝑡^𝑘_𝑖−𝑥^𝑡_𝑖+𝑐₂𝑟₂𝑔𝑏𝑒𝑠𝑡^𝑘_𝑖−𝑥^𝑘_𝑖,𝑥_𝑖^𝑘+1=𝑥^𝑘_𝑖+𝑣_𝑖^𝑘+1, (19) where ℎ represents the constriction factor and is defined as 2ℎ=2−𝛼−√𝛼²,−4𝛼𝛼=𝑐₁+𝑐₂>4. (20)(21)

In this variant of the PSO algorithm, ℎ controls the magnitude of the particle velocity and can be seen as a dampening factor. It provides the algorithm with two important features. First, it usually leads to faster convergence than standard PSO. Second, the swarm maintains the ability to perform wide movements in the search space, even if convergence is already advanced, but a new optimum is found. Therefore, the constriction PSO has the potential to avoid being trapped in local optima while possessing a fast convergence. It was shown to have superior performance compared to a standard PSO.

It is shown that a larger inertia weight tends to facilitate the global exploration and a smaller inertia weight achieves the local exploration to fine-tune the current search area [18]. The best performance could be obtained by initially setting 𝜔 to some relatively high value (e.g., 0.9), which corresponds to a system where particles perform extensive exploration, and gradually reducing 𝜔 to a much lower value (e.g., 0.4), where the system would be more dissipative and exploitative and would be better at homing into local optima. In [19], a linearly decreased inertia weight 𝜔 over time is proposed, where 𝜔 is given by the following equation: 𝜔=𝜔_𝑖−𝜔_𝑓⋅𝑘_max−𝑘𝑘_max+𝜔_𝑓, (22) where 𝜔_𝑖,𝜔_𝑓 are starting and final values of inertia weight, respectively; 𝑘_max is the maximum number of the iteration, and 𝑘 is the current iteration number. It is generally taken that starting value 𝜔_𝑖=0.9 and final value 𝜔_𝑓=0.4 [20].

On the other hand, in [21] was introduced PSO with time-varying acceleration coefficients. The improvement has the same motivation and the similar techniques as the adaptation of inertia weight. In this case, the cognitive coefficient 𝑐₁ is decreased linearly and the social coefficient 𝑐₂ is increased linearly over time as follows: 𝑐₁=𝑐_1𝑓−𝑐_1𝑖⋅𝑘_max−𝑘𝑘_max+𝑐_1𝑖,𝑐₂=𝑐_2𝑓−𝑐_2𝑖⋅𝑘_max−𝑘𝑘_max+𝑐_2𝑖, (23) where 𝑐_1𝑖 and 𝑐_2𝑖 are the initial values of the acceleration coefficients 𝑐₁ and 𝑐₂; 𝑐_1𝑓 and 𝑐_2𝑓 are the final values of the acceleration coefficients 𝑐₁ and 𝑐₂, respectively. Usually, 𝑐_1𝑖=2.5, 𝑐_2𝑖=0.5, 𝑐_1𝑓=0.5, and 𝑐_2𝑓=2.5.

3.3. Identification Algorithm

Considering all the states of the nonlinear system (8) at the sampling moments 𝑘∗𝑇𝑠  (𝑇𝑠 = sampling period) known, the identification algorithm has the following steps.

Step 1. Initialize a population of particles with random positions and velocities on 𝐷 dimensions in search space.

Step 2. For each particle, evaluate the desired optimization fitness function (17) in 𝐷 variables.

Step 3. Compare particle’s fitness evaluation with its pbest. If current value is better than pbest, then set pbest equal to the current value 𝑥_𝑖 in 𝐷-dimensional space.

Step 4. Identify the particle in swarm with the best success so far, and assign its index to the variable gbest.

Step 5. Change the velocity and position of the particle according to (19).

Step 6. if a criterion is met (usually a sufficiently good fitness or a maximum number of iterations)
then Stop;
else
go to Step 2.

3.4. “Classical” Identification Procedure

The “classical” approach for identification of yield coefficients is a two-step procedure under the assumption that full state measurements are available [1]. This method is based on a state transformation that allows reformulating the dynamical model into separate submodels. The first submodel depends only on the reaction structure and is independent of the kinetics. It can be linearly reparametrized and used for the identification of the yield coefficients by means of linear regressions, provided suitable identifiability conditions are satisfied. We present briefly this method for estimating the yield coefficients, and we use it for comparison in the next section.

The general dynamical model given in (2) represents a particular class of nonlinear state-space models. The nonlinearity lies in the reaction rates 𝜑_𝑖(𝜉) that are (nonlinear) functions of the state variables. These functions enter the model in the form 𝐾𝜑_𝑖(𝜉) (where 𝐾 is a constant matrix), which is a set of linear combinations of the same nonlinear functions 𝜑_𝑖(𝜉). This particular feature can be exploited to separate the nonlinear part from the linear part of the model by an adequate linear state transformation. More precisely, one chooses a nonsingular partition 𝐾_𝑎𝐾_𝑏=𝑇𝐾, (24) with 𝐾_𝑎∈𝑅^𝑝×𝑚 of full row rank matrix (i.e., 𝑝=rank(𝐾)), 𝐾_𝑏∈𝑅^{(𝑛−𝑝)×𝑚}, and 𝑇 a permutation matrix. The induced partitions of the vectors 𝜉 and 𝑢 are 𝜉_𝑎𝜉_𝑏𝑢=𝑇𝜉,^𝑎𝑢_𝑏=𝑇𝑢. (25) Model (2) is then partitioned into two submodels 𝑑𝜉_𝑎𝜉𝑑𝑡=𝐾⋅𝜑_𝑎,𝑡−𝐷𝜉_𝑎+𝑢_𝑎,𝑑𝜉_𝑏𝜉𝑑𝑡=𝐾⋅𝜑_𝑏,𝑡−𝐷𝜉_𝑏+𝑢_𝑏. (26) Then with the state transformation 𝜉_𝑎=𝜉_𝑎,𝑧=𝐶𝜉_𝑎+𝜉_𝑏, (27) one transforms the initial model into 𝜉̇_𝑎=𝐾_𝑎𝜑𝜉_𝑎,𝑧−𝐶𝜉_𝑎−𝐷𝜉_𝑎+𝑢_𝑎,𝑧̇=−𝐷𝑧+𝐶𝑢_𝑎+𝑢_𝑏, (28) where the (𝑛−𝑝)×𝑝 matrix 𝐶 is the unique solution of 𝐶𝐾_𝑎+𝐾_𝑏=0. (29) That is, 𝐶=−𝐾_𝑏𝐾^∗_𝑎, (30) where 𝐾^∗_𝑎 is a generalized inverse or pseudoinverse of 𝐾_𝑎. The subsystem (29) can be augmented with an equation derived from (8) as follows: 𝑧̇=−𝐷𝑧+𝐶𝑢_𝑎+𝑢_𝑏,𝜉_𝑏=𝑧−𝐶𝜉_𝑎. (31)

It can be considered as a linear time-varying (if 𝑑 varies in the course of time) model with state 𝑧, input (𝜉_𝑎, 𝑢_𝑎, and 𝑢_𝑏), and output 𝜉_𝑏. It is nonlinearly parametrized by the yield coefficients but linearly reparametrized by the nonzero entries of 𝐶. When data of the signals 𝜉_𝑎, 𝑢_𝑎, 𝑢_𝑏, and 𝜉_𝑏 are available, the auxiliary model (31) can be used to identify the yield coefficients independently of the knowledge of the reaction rates. The model (31) can be used to perform the identification of the nonzero entries of 𝐶 by a linear regression technique, with the yield coefficients 𝑘_𝑖 recovered afterwards from (29).

For the estimation of kinetic parameters the main approach is the use of a parameter observers (high gain observers, regressive parameter estimator, sliding mode observers, etc.). All these techniques have numerous tuning parameters and are difficult to implement.

4. Simulation Results and Discussion

The efficacy of our approach is shown by numerical simulations on an interval of 30 hours. The model given by relation (9) was integrated using a fourth-order Runge-Kutta routine with a sampling period of 1 minute and with initial conditions: IC=2.560.2100.50.

The influence of sampling period and type of the optimization algorithm and of noisy measurements are analyzed. To compare statistical performances of the different approaches the empirical normalized mean square error (NMSE) was used, that is defined as 1NMSE=𝑁^𝑁_𝑗=1NMSE𝜃_𝑗, (32) with NMSE(𝜃_𝑗)=((𝜃_𝑗−𝜃^∗_𝑗)/𝜃^∗_𝑗)², where 𝑁 is the number of estimated parameters, 𝜃_𝑗 is the 𝑗th element of the estimated parameter vector while the “*” superscript denotes the true value of the parameter.

In order to study the sensitivity of the estimation method to the sampling period and to the type of PSO algorithm, and to the noise, the following parameters were used.

Sampling period: 𝑇_𝑠∈{1min,10min,30min}.

Types of the optimization algorithm are as follows: PSO_1: algorithm based on relation (19) with ℎ defined by relation (20), and 𝑐₁=𝑐₂=2.1; PSO_2: algorithm based on relation (18) with 𝜔,𝑐₁, and 𝑐₂ defined by relations (22) and (23) and 𝜔_𝑖=0.9,𝜔_𝑓=0.4,𝑐_1𝑖=2.5;𝑐_2𝑖=0.5;𝑐_1𝑓=0.5,𝑐_2𝑓=2.5; PSO_3: algorithm based on relation (18) with 𝜔,𝑐₁,and𝑐₂ constants: 𝜔=0.75,𝑐₁=1.8,𝑐₂=2.2,  and 𝑉_max=10. Noise:{zero-mean white: SNR = 50 dB, 40 dB, 30 dB}.

The results of these simulations are presented in Tables 1, 2, and 3. The simulations were performed with a number of particles between 50 and 120. All the presented results are obtained for a number of 80 particles. For a greater number of particles the accuracy of the estimates was not better. In Figure 2 is presented the convergence rate for PSO_1 algorithm with a sampling period of 1 min.

Table 1 

Influence of the sampling period.

Table 2 

Influence of the algorithm type.

Table 3 

Influence of the noise level.

Figure 2 

Convergence rate for PSO_1 algorithm.

The results presented in Tables 1–3 suggest that our proposed parameter estimation technique yields consistent results. The method is very simple, easily completed and needs fewer parameters, which made it fully developed. However, the research on the PSO is still at the beginning, and a lot of problems are to be resolved. Research on the topology of the new pattern particle swarm which has a better function can be carried out. The neighbouring topology of the different particle swarms is based on the imitation of the different societies. It is meaningful to the use and spread of the algorithm to select the proper topology to enable PSO have the best property and do the research on the suitable ranges of different topologies. Blending PSO with the other intelligent optimization algorithms means combining the advantages of the PSO with the advantages of the other intelligent optimization algorithms to create the compound algorithm that has practical value. For example, the particle swarm optimization algorithm can be improved by the Simulated Annealing approach. Rapid swarm convergence is one of the main advantages of PSO, but this can also be problematic since, if an early solution is suboptimal, the swarm can easily stagnate around it without any pressure to continue exploration. Overall the results indicate that PSO algorithms can be used in the optimization of parameters during model identification.

5. Conclusions

The paper presents a particle swarm optimization based identification procedure for offline estimation of yield and kinetic coefficients in an anaerobic wastewater treatment bioprocess. The identification scheme is formulated in terms of an optimization problem where the error between an actual physical measured response of the system and the simulated response of a parameterized model is minimized. This function is multimodal, and classical iterative methods fail to find the global optimum. The estimation of the system parameters is achieved as a result of minimizing the error function by the PSO algorithm.

The simulations evaluated the effects of sampling period and some basic variants of PSO algorithm and of the noisy measurements. The proposed strategy can still converge to accurate results even in the presence of measurement noise, as illustrated by the numerical study. The PSO algorithm has a simpler procedure and higher computational efficiency than other optimization techniques.

Nomenclature

𝑐1,𝑐2:	Acceleration coefficients
𝐷:	Dilution rate (h−1)
𝐹:	Vector of feeding rates
𝐹in:	Input feed rate (Lh−1)
𝐾:	Matrix of yield coefficients
𝑘𝑖,𝑖=1,7:	Yield coefficients
𝑄:	Vector of rates of removal of the components in gaseous form
𝑆1:	Glucose concentration (g/L)
𝑆2:	Acetate concentration (g/L)
𝑆5:	Inorganic carbon concentration (g/L)
𝑋1:	Acidogenic concentration (g/L)
𝑋2:	Acetoclastic methanogenic concentration (g/L)
𝑃:	Methane concentration (g/L)
𝑆in:	Glucose concentration on the feed (g/L)
Φ:	Vector of reaction rates (reaction kinetics)
𝜇∗1,𝜇∗2:	Maximal specific growth rates (h−1)
𝜉:	State vector
ℜ:	The set of real numbers
Ω:	Inertia weight
𝜃:	The vector of unknown parameters
𝜃𝑖,𝑖=1,9:	The unknown parameters.

Acknowledgments

This work was supported by CNCSIS-UEFISCSU, Project no. PN II-RU TE 106, 9/04.08.2010.