# Parameter Identification of Anaerobic Wastewater Treatment Bioprocesses Using Particle Swarm Optimization

Authors: Dorin Sendrescu
Journal: Mathematical Problems in Engineering
 (2013)

---

## 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.


---

## Body

## 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 ofm biochemical reactions involving n components (with n≥m). The global dynamics can be represented by the following dynamical state-space model [1]:
(1)dξdt=K·φ(ξ,t)-Dξ+F-Q,
where ξ∈ℜn×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 K·φ(ξ,t) is the rate of consumption and/or production of the components in the reactor; that is, the reaction kinetics and the term -Dξ+F-Q 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 CO2, methane CH4, 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 H2, 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, H2, and inorganic carbon. In the first methanation path, acetate is transformed into CH4 and CO2 by acetoclastic methanogenic bacteria, while, in the second methanation path, H2 combines inorganic carbon to produce CH4 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:(2)X˙1=φ1-DX1,(3)S˙1=-k1φ1-DS1+DSin,(4)X˙2=φ2-DX2,(5)S˙2=k3φ1-k2φ2-DS2,(6)S˙5=k4φ1+k5φ2-DS5-QCO2,(7)P˙=k6φ2+k7φ1-DP-QP,
where X1, X2 are acidogenic and acetoclastic methanogenic bacteria, respectively, and S1, S2, S5 are glucose, acetate and inorganic carbon, respectively, and P is methane; φ1, φ2 are the rates of first acidogenic reaction and methanation reaction respectively, QP=cpP with cp>0 and QCO2 represent gaseous outflow rates of CH4 and CO2, respectively, Sin is the influent substrate concentration, D the dilution rate, and ki(i=0,…,7) are yield coefficients. Each reaction rate is a growth rate and may be written as φi=μiXi,i=1,2, where μi,i=1,2, is the specific growth rate of reaction i. Defining the state vector as ξ=[X1S1X2S2S5P]T, the model (2)–(7) can be written in matrix form as
(8)[X˙1S˙1X˙2S˙2S˙5P˙]=[10-k1001k2-k3k4k6k5k7][φ1φ2]-D[X1S1X2S2S5P]+[0DSin0000]-[0000QCO2QP],
or in compact form as
(9)ξ˙=Kφ(ξ)-Dξ+F-Q,
where F=[0DSin0000]T is the vector of inflow rates, Q=[0000QCO2QP]T is the vector of gaseous outflow rates, φ=[φ1φ2]T is the vector of reaction rates, which can be written as φ(ξ)=G(ξ)α(ξ), with G(ξ) being a diagonal matrix whose entries are products of the component concentrations involved in each reaction and α=[α1α2]T the vector of specific reaction rates, and K is the yield coefficient’s matrix. The matrices K and G have the following structure:
(10)K=[1-k10k2k4k6001-k3k5k7]T,G=[X1S100X2S2].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
(11)φ1(ξ)=μ1*S1·X1KM1+S1,
and the Haldane kinetic model
(12)φ2(ξ)=μ2*S2·X2KM2+S2+S22/Ki,
where KM1,KM2 are Michaelis-Menten constants, μ1*,μ2* represent specific growth rate coefficients, and Ki is the inhibition constant.For simplicity, we will denote the unknown plant parameters by the vector(13)θ=[θ1θ2⋯θ12]T,
where
(14)θ1=k1;θ2=k2;θ3=k3;θ4=k4;θ5=k5θ6=k6θ7=k7θ8=μ1*;θ9=μ2*;θ10=KM1;θ11=KM2;θ12=Ki.


## 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)dξ(t)dt=f(ξ,t;θ),
where ξ∈Rn is the state vector, θ∈Rm is the unknown parameters vector, and f  is a given nonlinear vector function.To estimate the unknown parameters in (15), a parameter identification system is defined as follows:
(16)dξ^(t)dt=f(ξ^,t;θ^),
where ξ^∈Rn is the estimated state vector and θ^∈Rm is the estimated parameters vector.The objective function defined as the mean squared errors between real and estimated responses for a numberN of given samples is considered as fitness of estimated model parameters:
(17)V=1N+D∑j=1D∑k=1N(ξjk-ξ^jk)2,
where D is the number of measurable states, N is the data length used for parameter identification, whereas ξjk and ξ^jk are the real and estimated values of state j at time k, 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 ofV 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 isD-dimensional, each individual is treated as a particle in the D-dimensional search space. The position and rate of position change for ith particle can be represented by D-dimensional vector, xi=(xi1,xi2,…,xiD) and vi=(vi1,vi2,…,viD), respectively. The best position previously visited by the ith particle is recorded and represented as pi=(pi1,pi2,…,piD), called pbest. The swarm best position previously visited by all the particles in the population is represented as pg=(pg1,pg2,…,pgD), called gbest. Then particles search their best position, which are guided by swarm information pg and their own information pi. 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:
(18)vik+1=ω·vik+c1r1(pbestik-xit)+c2r2(gbestik-xik),xik+1=xik+vik+1,
where k is iteration number, ω is inertia weight, c1 and c2 are two acceleration coefficients regulating the relative velocity toward local and global best positions, and r1 and r2  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 ω>(c1+c2)/2-1 and 0<ω<1. The technique originally proposed was to bound velocities so that each component of vi is kept within the range [-Vmax,+Vmax].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
(19)vik+1=h·[vik+c1r1(pbestik-xit)+c2r2(gbestik-xik)],xik+1=xik+vik+1,
where h represents the constriction factor and is defined as
(20)h=2|2-α-α2-4α|,(21)α=c1+c2>4.In this variant of the PSO algorithm,h 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:
(22)ω=(ωi-ωf)·kmax-kkmax+ωf,
where ωi,ωf are starting and final values of inertia weight, respectively; kmax is the maximum number of the iteration, and k is the current iteration number. It is generally taken that starting value ωi=0.9 and final value ωf=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 c1 is decreased linearly and the social coefficient  c2  is increased linearly over time as follows:
(23)c1=(c1f-c1i)·kmax-kkmax+c1i,c2=(c2f-c2i)·kmax-kkmax+c2i,
where c1i and c2i are the initial values of the acceleration coefficients c1 and c2; c1f and c2f are the final values of the acceleration coefficients c1 and c2, respectively. Usually, c1i=2.5, c2i=0.5, c1f=0.5, and c2f=2.5.


### 3.3. Identification Algorithm

Considering all the states of the nonlinear system (8) at the sampling moments k*Ts  (Ts = sampling period) known, the identification algorithm has the following steps.Step 1.
Initialize a population of particles with random positions and velocities onD dimensions in search space.Step 2.
For each particle, evaluate the desired optimization fitness function (17) in D variables.Step 3.
Compare particle’s fitness evaluation with itspbest. If current value is better than pbest, then set pbest equal to the current value xi in D-dimensional space.Step 4.
Identify the particle in swarm with the best success so far, and assign its index to the variablegbest.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 Step2.


### 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 φi(ξ) that are (nonlinear) functions of the state variables. These functions enter the model in the form Kφi(ξ) (where K is a constant matrix), which is a set of linear combinations of the same nonlinear functions φi(ξ). 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)(KaKb)=TK,
with Ka∈Rp×m of full row rank matrix (i.e., p=rank(K)), Kb∈R(n-p)×m, and T a permutation matrix. The induced partitions of the vectors ξ and u are
(25)(ξaξb)=Tξ,(uaub)=Tu.
Model (2) is then partitioned into two submodels
(26)dξadt=K·φ(ξa,t)-Dξa+ua,dξbdt=K·φ(ξb,t)-Dξb+ub.
Then with the state transformation
(27)ξa=ξa,z=Cξa+ξb,
one transforms the initial model into
(28)ξ˙a=Kaφ(ξa,z-Cξa)-Dξa+ua,z˙=-Dz+Cua+ub,
where the (n-p)×p matrix C is the unique solution of
(29)CKa+Kb=0.
That is,
(30)C=-KbKa*,
where Ka* is a generalized inverse or pseudoinverse of Ka. The subsystem (29) can be augmented with an equation derived from (8) as follows:
(31)z˙=-Dz+Cua+ub,ξb=z-Cξa.It can be considered as a linear time-varying (ifd varies in the course of time) model with state z, input (ξa, ua, and ub), and output ξb. It is nonlinearly parametrized by the yield coefficients but linearly reparametrized by the nonzero entries of C. When data of the signals ξa, ua, ub, and ξb 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 C by a linear regression technique, with the yield coefficients ki 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.


## 3.1. Problem Statement

Consider the following n-dimensional nonlinear system:(15)dξ(t)dt=f(ξ,t;θ),
where ξ∈Rn is the state vector, θ∈Rm is the unknown parameters vector, and f  is a given nonlinear vector function.To estimate the unknown parameters in (15), a parameter identification system is defined as follows:
(16)dξ^(t)dt=f(ξ^,t;θ^),
where ξ^∈Rn is the estimated state vector and θ^∈Rm is the estimated parameters vector.The objective function defined as the mean squared errors between real and estimated responses for a numberN of given samples is considered as fitness of estimated model parameters:
(17)V=1N+D∑j=1D∑k=1N(ξjk-ξ^jk)2,
where D is the number of measurable states, N is the data length used for parameter identification, whereas ξjk and ξ^jk are the real and estimated values of state j at time k, 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 ofV 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 isD-dimensional, each individual is treated as a particle in the D-dimensional search space. The position and rate of position change for ith particle can be represented by D-dimensional vector, xi=(xi1,xi2,…,xiD) and vi=(vi1,vi2,…,viD), respectively. The best position previously visited by the ith particle is recorded and represented as pi=(pi1,pi2,…,piD), called pbest. The swarm best position previously visited by all the particles in the population is represented as pg=(pg1,pg2,…,pgD), called gbest. Then particles search their best position, which are guided by swarm information pg and their own information pi. 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:
(18)vik+1=ω·vik+c1r1(pbestik-xit)+c2r2(gbestik-xik),xik+1=xik+vik+1,
where k is iteration number, ω is inertia weight, c1 and c2 are two acceleration coefficients regulating the relative velocity toward local and global best positions, and r1 and r2  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 ω>(c1+c2)/2-1 and 0<ω<1. The technique originally proposed was to bound velocities so that each component of vi is kept within the range [-Vmax,+Vmax].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
(19)vik+1=h·[vik+c1r1(pbestik-xit)+c2r2(gbestik-xik)],xik+1=xik+vik+1,
where h represents the constriction factor and is defined as
(20)h=2|2-α-α2-4α|,(21)α=c1+c2>4.In this variant of the PSO algorithm,h 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:
(22)ω=(ωi-ωf)·kmax-kkmax+ωf,
where ωi,ωf are starting and final values of inertia weight, respectively; kmax is the maximum number of the iteration, and k is the current iteration number. It is generally taken that starting value ωi=0.9 and final value ωf=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 c1 is decreased linearly and the social coefficient  c2  is increased linearly over time as follows:
(23)c1=(c1f-c1i)·kmax-kkmax+c1i,c2=(c2f-c2i)·kmax-kkmax+c2i,
where c1i and c2i are the initial values of the acceleration coefficients c1 and c2; c1f and c2f are the final values of the acceleration coefficients c1 and c2, respectively. Usually, c1i=2.5, c2i=0.5, c1f=0.5, and c2f=2.5.


## 3.3. Identification Algorithm

Considering all the states of the nonlinear system (8) at the sampling moments k*Ts  (Ts = sampling period) known, the identification algorithm has the following steps.Step 1.
Initialize a population of particles with random positions and velocities onD dimensions in search space.Step 2.
For each particle, evaluate the desired optimization fitness function (17) in D variables.Step 3.
Compare particle’s fitness evaluation with itspbest. If current value is better than pbest, then set pbest equal to the current value xi in D-dimensional space.Step 4.
Identify the particle in swarm with the best success so far, and assign its index to the variablegbest.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 Step2.


## 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 φi(ξ) that are (nonlinear) functions of the state variables. These functions enter the model in the form Kφi(ξ) (where K is a constant matrix), which is a set of linear combinations of the same nonlinear functions φi(ξ). 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)(KaKb)=TK,
with Ka∈Rp×m of full row rank matrix (i.e., p=rank(K)), Kb∈R(n-p)×m, and T a permutation matrix. The induced partitions of the vectors ξ and u are
(25)(ξaξb)=Tξ,(uaub)=Tu.
Model (2) is then partitioned into two submodels
(26)dξadt=K·φ(ξa,t)-Dξa+ua,dξbdt=K·φ(ξb,t)-Dξb+ub.
Then with the state transformation
(27)ξa=ξa,z=Cξa+ξb,
one transforms the initial model into
(28)ξ˙a=Kaφ(ξa,z-Cξa)-Dξa+ua,z˙=-Dz+Cua+ub,
where the (n-p)×p matrix C is the unique solution of
(29)CKa+Kb=0.
That is,
(30)C=-KbKa*,
where Ka* is a generalized inverse or pseudoinverse of Ka. The subsystem (29) can be augmented with an equation derived from (8) as follows:
(31)z˙=-Dz+Cua+ub,ξb=z-Cξa.It can be considered as a linear time-varying (ifd varies in the course of time) model with state z, input (ξa, ua, and ub), and output ξb. It is nonlinearly parametrized by the yield coefficients but linearly reparametrized by the nonzero entries of C. When data of the signals ξa, ua, ub, and ξb 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 C by a linear regression technique, with the yield coefficients ki 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(32)NMSE=1N∑j=1NNMSE(θ^j),
with NMSE(θ^j)=((θ^j-θj*)/θj*)2, where N is the number of estimated parameters, θ^j is the jth 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:Ts∈{1min,10min,30min}.Types of the optimization algorithm are as follows:PSO_1:
algorithm based on relation (19) with h defined by relation (20), and c1=c2=2.1;PSO_2:
algorithm based on relation (18) with ω,c1, and c2 defined by relations (22) and (23) and ωi=0.9,ωf=0.4,c1i=2.5;c2i=0.5;c1f=0.5,c2f=2.5;PSO_3:
algorithm based on relation (18) with ω,c1,andc2 constants: ω=0.75,c1=1.8,c2=2.2,  and Vmax=10.Noise:
{zero-mean white: SNR = 50 dB, 40 dB,  30 dB}.The results of these simulations are presented in Tables1, 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.
Estimated parameter
“Real” value
Sampling period
1 min
10 min
30 min
θ
1
5.5
5.5009
5.4811
5.4798
θ
2
1
0.9982
1.0111
1.0637
θ
3
16
15.9110
15.9328
14.5629
θ
4
10
10.0016
9.9707
9.9571
θ
5
1.5
1.4938
1.4626
1.3078
θ
6
0.2
0.2003
0.2011
0.1873
θ
7
3
2.9823
2.9581
2.6947
θ
8
0.2
0.2002
0.1974
0.1949
θ
9
0.6
0.5994
0.6010
0.6188
θ
10
0.75
0.7539
0.6754
0.6483
θ
11
4
3.9731
3.9007
3.0265
θ
12
20
20.3973
19.4738
20.3947
NMSE
0
3.0717e − 05
5.6623e − 04
1.2143e − 01Table 2
Influence of the algorithm type.
Estimated parameter
“Real” value
Type of optimization algorithm
PSO-1
PSO-2
PSO-3
θ
1
5.5
5.5613
5.5849
5.4338
θ
2
1
1.0405
1.0083
0.8406
θ
3
16
15.2050
14.9971
15.6475
θ
4
10
10.0848
10.1252
9.9468
θ
5
1.5
1.5721
1.6422
1.1698
θ
6
0.2
0.2022
0.2040
0.3005
θ
7
3
2.8006
2.7916
2.7569
θ
8
0.2
0.1990
0.1975
0.2286
θ
9
0.6
0.5993
0.6311
0.6594
θ
10
0.75
0.9991
1.1153
0.4728
θ
11
4
3.5221
3.3470
4.0707
θ
12
20
20.9257
20.5048
20.0971
NMSE
0
3.0717e − 04
0.0283
0.3845Table 3
Influence of the noise level.
Estimated parameter
“Real” value
Noise level
50 dB
40 dB
30 dB
θ
1
5.5
5.4001
5.3697
5.3362
θ
2
1
1.0066
0.6628
0.4180
θ
3
16
15.5721
16.6508
16.3867
θ
4
10
9.8502
9.8463
9.7771
θ
5
1.5
1.1817
1.1428
0.9150
θ
6
0.2
0.1867
0.2769
0.2711
θ
7
3
2.9626
2.9178
3.1112
θ
8
0.2
0.2042
0.2350
0.3020
θ
9
0.6
0.6030
0.6254
0.5744
θ
10
0.75
0.3306
0.2379
0.0058
θ
11
4
3.7860
4.7660
4.5629
θ
12
20
20.8802
20.9267
19.7113
NMSE
0
3.6827e − 3
0.0827
0.1625Figure 2
Convergence rate for PSO_1 algorithm.The results presented in Tables1–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.



---

*Source: 103748-2013-07-17.xml*