# Parameter Identification of the 2-Chlorophenol Oxidation Model Using Improved Differential Search Algorithm

Authors: Guang-zhou Chen; Jia-quan Wang; Ru-zhong Li
Journal: Journal of Chemistry
 (2015)

---

## Abstract

Parameter identification plays a crucial role for simulating and using model. This paper firstly carried out the sensitivity analysis of the 2-chlorophenol oxidation model in supercritical water using the Monte Carlo method. Then, to address the nonlinearity of the model, two improved differential search (DS) algorithms were proposed to carry out the parameter identification of the model. One strategy is to adopt the Latin hypercube sampling method to replace the uniform distribution of initial population; the other is to combine DS with simplex method. The results of sensitivity analysis reveal the sensitivity and the degree of difficulty identified for every model parameter. Furthermore, the posteriori probability distribution of parameters and the collaborative relationship between any two parameters can be obtained. To verify the effectiveness of the improved algorithms, the optimization performance of improved DS in kinetic parameter estimation is studied and compared with that of the basic DS algorithm, differential evolution, artificial bee colony optimization, and quantum-behaved particle swarm optimization. And the experimental results demonstrate that the DS with the Latin hypercube sampling method does not present better performance, while the hybrid methods have the advantages of strong global search ability and local search ability and are more effective than the other algorithms.


---

## Body

## 1. Introduction

In environmental technology, supercritical water oxidation (SCWO) is an innovative technology. Supercritical water exhibits some properties very different from those of water at standard conditions. One of its main applications is the oxidative treatment of aqueous wastes containing organic compounds in the so-called “supercritical water oxidation” (SCWO) process [1]. Its potential advantage lies in the fact that it can rapidly and efficiently decontaminate diverse organic wastes [2]. In wastewater treatment, supercritical water oxidation technology is considered to be an effective clean technology and is applied to the organic wastewater treatment [3–5].The key technical issue in the design and application of SCWO process systems is to understand the principle of kinetic reaction and identify the key reaction mechanisms and parameters in the kinetic model [2]. In modeling, the mathematical model is designed to simulate the system behavior f
(
x
) for a given set of model parameters x representing the actual response properties of the system. The more accurately the model matches the experiment observations, the more representative the model parameters are assumed to be. Therefore, the above parameter identification or parameter estimation is essentially an inverse problem [6].In many applications, inverse problems can be nonlinear and ill posed which makes them difficult to solve numerically. The basic idea is to carry out the optimization by constructing the objective function. The best set of parameters is selected from within reasonable ranges by adjusting the values until the discrepancies between the model generated values and observations are reduced to a minimum in the weighted least squares sense [7]. In addition, in the process of model identification, along with calibration to obtain a set of optimal parameter sets giving the best performance of a simulation model, there are other important issues such as uncertainty and identify ability of the parameters which should be taken into account for a proper model evaluation.Generally speaking, the optimization methods for solving the inverse problems are classified into two groups. One sort is traditional methods: (1) gradient-based optimal methods such as Levenberg-Marquardt algorithm and the conjugate gradient descent method; (2) no gradient-based algorithms such as simplex. And the other is search-based heuristic optimization methods such as genetic algorithm (GA), simulated annealing (SA), and tabu search algorithm (TA), which are used in dealing with complex nonlinear inverse problems. They have been proved to have the ability to identify global or near-optimal solutions in the global scope by a large number of engineering optimization cases.In recent years, many swarm intelligence optimization algorithms, such as ant colony, particle swarm optimization, artificial fish swarm, and artificial bee colony, have been put forward unceasingly. Many swarm intelligence algorithms have been applied in the parameter identification from diverse engineering domains. For the chemical kinetic models, chaos genetic algorithm, clonal selection algorithm, hierarchical differential evolution algorithm, adaptive differential evolution algorithm, and ABC-simplex and DNA genetic algorithm have been applied to estimate the kinetic parameters [8–13].Though the above intelligence optimization algorithms can usually obtain better optimization results, a more effective and convenient optimization algorithm is still worth seeking. Differential search (DS) algorithm, proposed by Civicioglu in 2012 [14], is a new optimization algorithm with a short code, short run time, and simple principle. Subsequently, some applications show that it is an effective optimization algorithm and has better optimization performance [14, 15].Aimed at the case where the sensitive and uncertainty analysis is easy to be ignored in parameter identification process, the main objectives of this paper are to use Monte Carlo method to carry out sensitivity analysis based on Monte Carlo analysis toolbox (MCAT); to develop two improved algorithms to enhance optimization performance; to analyze their optimization performance by a case and compare it with that of some other swarm intelligence optimization algorithms.This paper is organized as follows. In Section2, we introduce the sensitivity analysis and MCAT. Then, the improved DS algorithms are presented in Section 3. In Section 4, we apply MCAT to make sensitivity analysis and improved DS algorithms to test a case of parameter identification from a chemical kinetics model and analyze its performance. Finally, we give some conclusions in Section 5.


## 2. Sensitive Analysis and Introduction of MCAT

### 2.1. Sensitive Analysis

To provide guidance for a better understanding of the modeling and reduce the response uncertainties most effectively, sensitivity analysis (SA) of the input importance on the response variability can be useful [16]. Sensitivity analysis is valuable tool for identifying important model parameters, testing the model conceptualization, and improving the model structure. They help to apply the model efficiently and to enable a focused planning of future research and field measurement [17].


### 2.2. Monte Carlo Analysis Toolbox

The Monte Carlo analysis toolbox (MCAT), developed by Wagener et al., is a collection of analysis and visualization functions integrated through a graphical user interface based on MATLAB language [18]. The toolbox can be used to analyze the results from Monte Carlo (MC) parameter sampling experiments or from model optimization methods that are based on population evolution techniques. A number of powerful functions are included in the toolbox to investigate the structure, sensitivity, and parameter and output uncertainty of models. In view of the problem that system identification techniques are often limited to linear systems and are subject to many statistical assumptions, MCAT uses the Monte Carlo simulation approach to solve the “equifinality” [19] instead of classical system identification techniques [20]. Monte Carlo parameter sampling provides an alternative measure to classical statistical modeling and uncertainty estimation methods. It is designed to help interpret the results of such modeling experiments.The core of the toolbox is based on the concept of Regional Sensitivity Analysis [21–23], and its extension to the generalized likelihood uncertainty estimation (GLUE) technique was developed at Lancaster University by Beven and coworkers [19, 23]. For detailed description, see guide of MCAT.


## 2.1. Sensitive Analysis

To provide guidance for a better understanding of the modeling and reduce the response uncertainties most effectively, sensitivity analysis (SA) of the input importance on the response variability can be useful [16]. Sensitivity analysis is valuable tool for identifying important model parameters, testing the model conceptualization, and improving the model structure. They help to apply the model efficiently and to enable a focused planning of future research and field measurement [17].


## 2.2. Monte Carlo Analysis Toolbox

The Monte Carlo analysis toolbox (MCAT), developed by Wagener et al., is a collection of analysis and visualization functions integrated through a graphical user interface based on MATLAB language [18]. The toolbox can be used to analyze the results from Monte Carlo (MC) parameter sampling experiments or from model optimization methods that are based on population evolution techniques. A number of powerful functions are included in the toolbox to investigate the structure, sensitivity, and parameter and output uncertainty of models. In view of the problem that system identification techniques are often limited to linear systems and are subject to many statistical assumptions, MCAT uses the Monte Carlo simulation approach to solve the “equifinality” [19] instead of classical system identification techniques [20]. Monte Carlo parameter sampling provides an alternative measure to classical statistical modeling and uncertainty estimation methods. It is designed to help interpret the results of such modeling experiments.The core of the toolbox is based on the concept of Regional Sensitivity Analysis [21–23], and its extension to the generalized likelihood uncertainty estimation (GLUE) technique was developed at Lancaster University by Beven and coworkers [19, 23]. For detailed description, see guide of MCAT.


## 3. Improved Differential Search Algorithm

### 3.1. Main Idea of the Basic DS Algorithm

DS is an algorithm developed for solution of numeric optimization problems [14]. It was inspired by migration of living beings which constitute superorganisms utilizing the concept of Brownian-like motion. Suffering from the impact of periodical change of regional resources, many species of the living beings show seasonal migration behavior during the year. In the migration movement, the migrating species of living beings constitute a superorganism containing large number of individuals.Migration behavior allows them to move from one habitat to more efficient habitat. Of course, they change their position by moving toward more fruitful areas. The movement of superorganisms can be described by a Brownian-like rand-walk model [24]. In algorithm, evolving process of populations is corresponding to artificial superorganisms migration to global optimum solution of the problem. During the migration, the superorganisms judge whether some randomly selected positions are suitable for temporary stopovers. If the position is suitable to stop over for a temporary time during the migration, the members of superorganisms that made the discovery settle at this position; subsequently, they continue their migration from this position.


### 3.2. Main Steps of the Basic DS

Step 1. Set up the size of population, dimension of the problem, and numbers of iterations.
Step 2. Evaluate the initial population.
Step 3. Execution cycle is as follows:
(1)
use the random shuffling method to produce the donor population;
(2)
compute the scale factor;
(3)
obtain the stopover site position population;
(4)
use the random process to decide the individuals to participate in the search process of stopover site;
(5)
check whether the elements of stopover site are beyond the limits of the search space; if beyond, a position in the search space will be randomly produced;
(6)
evaluate the individuals of the stopover site position population; if it is better, replace the original one; or it remains unchanged;
(7)
update the entire group.
Step 4. If the numbers of iterations are satisfied, then stop; otherwise, continue step 3.


### 3.3. Improved DS Algorithms

#### 3.3.1. Construction of the Initial Population Using Latin Hypercube Sampling Method

Latin hypercube sampling (LHS) is a statistical method for generating a sample of plausible collections of parameter values from a multidimensional distribution. The sampling method is often used to construct the design of experiments. Latin hypercube sampling can be viewed as a compromise procedure that incorporates many of the desirable features of random sampling and stratified sampling and also produces more stable analysis outcomes than random sampling [25]. This sampling approach ensures that each of the input variables has all portions of its range represented; LHS is computationally cheap to generate and can cope with many input variables [26]. Consequently, LHS was used to produce the initial population instead of uniform distribution method so as to enhance the performance of algorithm. Its realization was simply described [26]: generate a sample size N from the n variables ξ
1
,
ξ
2
,
…
,
ξ
n. The range of each variable is divided into N nonoverlapping intervals on the basis of equal probability size 1
/
N. One value from each interval is selected at random with respect to the probability density in the interval. The N values thus obtained for ξ
1 are paired in a random manner with the N values of ξ
2. These N pairs are combined randomly with the N values of ξ
3 to form N triplets and so on, until a set of N  n-tuples is formed. This set of n-tuples is the Latin hypercube sample. In this paper, we use “lhsdesign” function of MATLAB to realize it.


#### 3.3.2. Combination of DS and Simplex

Simplex was proposed by Spendley et al. in 1962. Subsequently, it was improved by Nelder and Mead in 1965 [27]. It is a direct search optimization method, and has the fast searching ability to the local optimization without the need of gradient information. Its biggest weakness is sensitive to the choice of the initial points. For the DS algorithm, though it has the powerful global search ability, the convergence speed of the basic DS algorithm is slow and needs a large number of iterations.Therefore, to own both the merits of the better local searching ability from the simplex method and the global searching ability from the DS algorithm, the hybrid algorithms of DS and simplex are developed.The paper designed two methods: (1) firstly completely run the DS algorithm to obtain the optimization solution; then simplex method is employed to process local search based on the above solution which is considered as the initial solution; (2) at the end of each run of the DS algorithm, simplex method is employed to execute local search based on the optimal solution of the current generation. If the better solution can be obtained, the optimal solution of the current generation will be replaced. Therefore, it effectively absorbs the advantages of DS and simplex method. It not only can overcome the deficiency that DS has weak local search ability but also can effectively solve the sensitivity problem of simplex to the initial search points.


## 3.1. Main Idea of the Basic DS Algorithm

DS is an algorithm developed for solution of numeric optimization problems [14]. It was inspired by migration of living beings which constitute superorganisms utilizing the concept of Brownian-like motion. Suffering from the impact of periodical change of regional resources, many species of the living beings show seasonal migration behavior during the year. In the migration movement, the migrating species of living beings constitute a superorganism containing large number of individuals.Migration behavior allows them to move from one habitat to more efficient habitat. Of course, they change their position by moving toward more fruitful areas. The movement of superorganisms can be described by a Brownian-like rand-walk model [24]. In algorithm, evolving process of populations is corresponding to artificial superorganisms migration to global optimum solution of the problem. During the migration, the superorganisms judge whether some randomly selected positions are suitable for temporary stopovers. If the position is suitable to stop over for a temporary time during the migration, the members of superorganisms that made the discovery settle at this position; subsequently, they continue their migration from this position.


## 3.2. Main Steps of the Basic DS

Step 1. Set up the size of population, dimension of the problem, and numbers of iterations.
Step 2. Evaluate the initial population.
Step 3. Execution cycle is as follows:
(1)
use the random shuffling method to produce the donor population;
(2)
compute the scale factor;
(3)
obtain the stopover site position population;
(4)
use the random process to decide the individuals to participate in the search process of stopover site;
(5)
check whether the elements of stopover site are beyond the limits of the search space; if beyond, a position in the search space will be randomly produced;
(6)
evaluate the individuals of the stopover site position population; if it is better, replace the original one; or it remains unchanged;
(7)
update the entire group.
Step 4. If the numbers of iterations are satisfied, then stop; otherwise, continue step 3.


## 3.3. Improved DS Algorithms

### 3.3.1. Construction of the Initial Population Using Latin Hypercube Sampling Method

Latin hypercube sampling (LHS) is a statistical method for generating a sample of plausible collections of parameter values from a multidimensional distribution. The sampling method is often used to construct the design of experiments. Latin hypercube sampling can be viewed as a compromise procedure that incorporates many of the desirable features of random sampling and stratified sampling and also produces more stable analysis outcomes than random sampling [25]. This sampling approach ensures that each of the input variables has all portions of its range represented; LHS is computationally cheap to generate and can cope with many input variables [26]. Consequently, LHS was used to produce the initial population instead of uniform distribution method so as to enhance the performance of algorithm. Its realization was simply described [26]: generate a sample size N from the n variables ξ
1
,
ξ
2
,
…
,
ξ
n. The range of each variable is divided into N nonoverlapping intervals on the basis of equal probability size 1
/
N. One value from each interval is selected at random with respect to the probability density in the interval. The N values thus obtained for ξ
1 are paired in a random manner with the N values of ξ
2. These N pairs are combined randomly with the N values of ξ
3 to form N triplets and so on, until a set of N  n-tuples is formed. This set of n-tuples is the Latin hypercube sample. In this paper, we use “lhsdesign” function of MATLAB to realize it.


### 3.3.2. Combination of DS and Simplex

Simplex was proposed by Spendley et al. in 1962. Subsequently, it was improved by Nelder and Mead in 1965 [27]. It is a direct search optimization method, and has the fast searching ability to the local optimization without the need of gradient information. Its biggest weakness is sensitive to the choice of the initial points. For the DS algorithm, though it has the powerful global search ability, the convergence speed of the basic DS algorithm is slow and needs a large number of iterations.Therefore, to own both the merits of the better local searching ability from the simplex method and the global searching ability from the DS algorithm, the hybrid algorithms of DS and simplex are developed.The paper designed two methods: (1) firstly completely run the DS algorithm to obtain the optimization solution; then simplex method is employed to process local search based on the above solution which is considered as the initial solution; (2) at the end of each run of the DS algorithm, simplex method is employed to execute local search based on the optimal solution of the current generation. If the better solution can be obtained, the optimal solution of the current generation will be replaced. Therefore, it effectively absorbs the advantages of DS and simplex method. It not only can overcome the deficiency that DS has weak local search ability but also can effectively solve the sensitivity problem of simplex to the initial search points.


## 3.3.1. Construction of the Initial Population Using Latin Hypercube Sampling Method

Latin hypercube sampling (LHS) is a statistical method for generating a sample of plausible collections of parameter values from a multidimensional distribution. The sampling method is often used to construct the design of experiments. Latin hypercube sampling can be viewed as a compromise procedure that incorporates many of the desirable features of random sampling and stratified sampling and also produces more stable analysis outcomes than random sampling [25]. This sampling approach ensures that each of the input variables has all portions of its range represented; LHS is computationally cheap to generate and can cope with many input variables [26]. Consequently, LHS was used to produce the initial population instead of uniform distribution method so as to enhance the performance of algorithm. Its realization was simply described [26]: generate a sample size N from the n variables ξ
1
,
ξ
2
,
…
,
ξ
n. The range of each variable is divided into N nonoverlapping intervals on the basis of equal probability size 1
/
N. One value from each interval is selected at random with respect to the probability density in the interval. The N values thus obtained for ξ
1 are paired in a random manner with the N values of ξ
2. These N pairs are combined randomly with the N values of ξ
3 to form N triplets and so on, until a set of N  n-tuples is formed. This set of n-tuples is the Latin hypercube sample. In this paper, we use “lhsdesign” function of MATLAB to realize it.


## 3.3.2. Combination of DS and Simplex

Simplex was proposed by Spendley et al. in 1962. Subsequently, it was improved by Nelder and Mead in 1965 [27]. It is a direct search optimization method, and has the fast searching ability to the local optimization without the need of gradient information. Its biggest weakness is sensitive to the choice of the initial points. For the DS algorithm, though it has the powerful global search ability, the convergence speed of the basic DS algorithm is slow and needs a large number of iterations.Therefore, to own both the merits of the better local searching ability from the simplex method and the global searching ability from the DS algorithm, the hybrid algorithms of DS and simplex are developed.The paper designed two methods: (1) firstly completely run the DS algorithm to obtain the optimization solution; then simplex method is employed to process local search based on the above solution which is considered as the initial solution; (2) at the end of each run of the DS algorithm, simplex method is employed to execute local search based on the optimal solution of the current generation. If the better solution can be obtained, the optimal solution of the current generation will be replaced. Therefore, it effectively absorbs the advantages of DS and simplex method. It not only can overcome the deficiency that DS has weak local search ability but also can effectively solve the sensitivity problem of simplex to the initial search points.


## 4. Parameter Estimation for Kinetic Model of the 2-Chlorophenol

### 4.1. The Kinetic Model

Supercritical water oxidation (SCWO) is an effective treatment technology for organic waste in environmental field. The removing rate expression for 2-chlorophenol (2-CP) is as follows [28]:
(1)
r
=
A
exp
⁡
-
E
a
R
T
2
CP
a
O
2
b
H
2
O
c
,
where [
2
CP
], [
O
2
], and [
H
2
O
] represent 2-CP, O
2, and H
2
O concentration, respectively. r is the removing rate of 2-CP, A is the preexponential factor, E
a denotes the activation energy, R is the molar gas constant, and a, b, and c are the reaction orders of 2-CP,  O
2, and H
2
O, respectively. After the analysis and treatment, the following equations can be obtained [28]:
(2)
ln
⁡
1
-
X
=
-
A
exp
⁡
-
E
a
R
T
O
2
b
H
2
O
c
τ
if
a
=
1
,
(3)
1
-
X
1
-
a
-
1
=
a
-
1
A
exp
⁡
-
E
a
R
T
2
CP
a
-
1
O
2
b
H
2
O
c
τ
if
a
≠
1
,
where X represents the conversion rate of 2-CP. The objective is to identify the five parameters (A
,
E
a
,
a
,
b
,
c). This optimization problem is a complex one with high nonlinearity.


### 4.2. The Objective Function

This parameter estimation problem is based on minimizing the objective function given by the sum of squared deviation. The objective function is defined as(4)
min
⁡
⁡
f
=
∑
i
=
1
n
X
i
-
X
i
′
2
,
where n is the sample number and X
i is the real conversion rate of 2-CP in the ith sample data while X
i
′ is the conversion rate calculated by the kinetic model.


### 4.3. Sensitivity Analysis

We used MCAT to carry out the sensitivity analysis and uncertainty analysis. After the 50000 runs, the threshold of sum of squared errors (SSE) for the objective function was taken to 2, and then the following plots were obtained.


#### 4.3.1. Dotty Plots of Objective Function

The dotty plots shown in Figure1 is the result of rejection threshold of sum of squared errors (SSE) set to 2. The “surfaces” of the dotty plot for a, b, and c parameters have a clearly defined minimum, so they can be considered to be well identified. Followed by E
a, however, the A parameter is difficult to be defined.Figure 1
Dotty plots of objective function (SSE).


#### 4.3.2. Posteriori Parameter Distribution

The initial simulation population was uniformly distributed, and the resulting distribution valuesD are plotted as bars as shown in Figure 2. In plots, the range of every parameter is divided into 20 containers of equal width. The likelihoods in each container are added up and divided by the sum of parameter values within the container. The regions of higher bars in Figure 2 show that their corresponding parameter ranges have better model performance. So, we can obtain the better region of parameters aimed at the model from Figure 2.Figure 2
Posteriori parameter distribution.


#### 4.3.3. Regional Sensitivity Analysis Plot

The parameter sets are then split into 10 groups. For each group the likelihoods are normalized by dividing by their total, and the cumulative frequency distribution is calculated and plotted. If the model performance is sensitive to a particular parameter, there will be a large difference between the cumulative frequency distributions. Figure3 shows that parameter c is most sensitive, followed by parameters a, b, and E
a. And parameter A is insensitive.Figure 3
Regional sensitivity analysis plot.


#### 4.3.4. 3D Surface Plot of Relationship between the Two Parameters

In order to examine parameter interaction, a Monte Carlo simulation was run with just 2 parameters varying and the other parameters fixed for MCAT. The plots in Figure4 show the parameter surfaces of every two parameters. As seen in Figure 4, we can note that there is a relatively well-defined region of low SSE, which indicates better combination state of the parameters. Therefore, it can provide the basis to judge the results of parameter identification.Figure 4
3D surface plot of every two parameters.
(a)
A
-
E
a
(b)
A-a
(c)
A-b
(d)
A-c
(e)
E
a-a
(f)
E
a-b
(g)
E
a-c
(h)
a-b
(i)
a-c
(j)
b-c


#### 4.3.5. GLUE Variable Uncertainty

Figure5 shows the cumulative probability distribution and probability density function for the output variable calculated using a selected objective (transformed to likelihood). It also shows the uncertainty in peak output for the above example. The result shows that there is a probability of 0.05 (5%) that the peak output will be above 0.999 or below 0.865.Figure 5
GLUE variable uncertainty plot.The above results can help us to identify the important model input parameters, identify the main factors of uncertainty for model output, improve the model structure, and judge the rationality of identification results.


### 4.4. Parameter Estimation

62 groups of experimental data including temperatureT, pressure P, residence time τ, [
2
CP
], [
O
2
], [
H
2
O
], and 2-CP conversion X provided in the literature [28] are used as test samples to verify the efficiency of the proposed algorithms.In the experiments, the population of DS is set to 50, the total iterative number is 500, scale factor is produced by using the gamma distribution, others are set with default parameters, and the parameter ranges of are[
0,200
], [
30000,50000
], [
0,2
], [
0,2
], and [
0,2
], respectively. Taking into account the situation in which DS has shorter running time, we increased an experiment where the total iterative number is 10000.As a comparison, the population of QPSO is set to 50; the total iterative number is 500;α decreases linearly from 1 to 0.5; the population of bees in ABC is set to 50; the total iterative number is 500; the limit value is 100; the ranges of parameters are the same as the previous ones.In order to prevent influence of random factor, the above algorithms are executed in 20 independent runs and the best, average, and the worst values of the results are recorded. The optimization results of the above algorithms are listed in Table1 in detail. As a comparison, the best results from some other algorithms are as follows: 0.2225 (chaos GA), 0.2494 (nonlinear regression), 0.2177 (adaptive differential evolution), 0.2177 (hybrid ant colony system), 0.2753 (RNA-GA, [13]), and 0.2551 (MDNA-GA, [13]).Table 1
Objective function values obtained by different algorithms.
Algorithms
Best
Average
Worst
DE [10]
0.217685
0.217725
0.218037
PSO [10]
0.217685
0.217819
0.218733
Hierarchical differential evolution [10]
0.217685
0.217685
0.217685
ABC [12]
0.233919
0.292802
0.421080
Hybrid ABC (HABC) [12]
0.217685
0.217685
0.217685
QPSO [29]
0.2176855
0.2177133
0.2178465
DS(1)
0.222266
0.252299
0.287777
DS(2)
0.2176850
0.2176855
0.217689
DS(lhs)
0.222360
0.248886
0.308399
Improved DS(1)
0.2176850
0.2176850
0.2176850
Improved DS(2)
0.2176850
0.2176850
0.2176850
Note: DS(1) represents the iteration number which is 500; DS(2) represents the iteration number which is 10000; DS(lhs) represents the initial population which is generated by Latin hypercube sampling; Improved DS(1) represents the simplex method which is used according to the sequential order; Improved DS(2) represents the other hybrid method, whose total iteration number is 5.Table1 shows that the best objective function values achieved by the algorithms (DE, PSO, HDE, HABC, QPSO, DS(2), Improved DS(1), and Improved DS(2)) are very similar, while the average and worst values obtained by HDE, hybrid ABC, DS(2), and the two improved DS algorithms outperform those of DE, PSO, and ABC. Furthermore, the basic DS algorithm can be used successfully with a large number of runs, and DS(lhs) does not show better optimization results. The obtained results also show that the two improvement algorithms compared with simplex obtain better performance, which has the advantages of strong global search ability and local search ability, can complete the optimization task outstandingly. The optimum solution (A
,
E
a
,
a
,
b
,
c) is 63.5468, 45626, 0.8081, 0.4444, and 0.3239, respectively.


## 4.1. The Kinetic Model

Supercritical water oxidation (SCWO) is an effective treatment technology for organic waste in environmental field. The removing rate expression for 2-chlorophenol (2-CP) is as follows [28]:
(1)
r
=
A
exp
⁡
-
E
a
R
T
2
CP
a
O
2
b
H
2
O
c
,
where [
2
CP
], [
O
2
], and [
H
2
O
] represent 2-CP, O
2, and H
2
O concentration, respectively. r is the removing rate of 2-CP, A is the preexponential factor, E
a denotes the activation energy, R is the molar gas constant, and a, b, and c are the reaction orders of 2-CP,  O
2, and H
2
O, respectively. After the analysis and treatment, the following equations can be obtained [28]:
(2)
ln
⁡
1
-
X
=
-
A
exp
⁡
-
E
a
R
T
O
2
b
H
2
O
c
τ
if
a
=
1
,
(3)
1
-
X
1
-
a
-
1
=
a
-
1
A
exp
⁡
-
E
a
R
T
2
CP
a
-
1
O
2
b
H
2
O
c
τ
if
a
≠
1
,
where X represents the conversion rate of 2-CP. The objective is to identify the five parameters (A
,
E
a
,
a
,
b
,
c). This optimization problem is a complex one with high nonlinearity.


## 4.2. The Objective Function

This parameter estimation problem is based on minimizing the objective function given by the sum of squared deviation. The objective function is defined as(4)
min
⁡
⁡
f
=
∑
i
=
1
n
X
i
-
X
i
′
2
,
where n is the sample number and X
i is the real conversion rate of 2-CP in the ith sample data while X
i
′ is the conversion rate calculated by the kinetic model.


## 4.3. Sensitivity Analysis

We used MCAT to carry out the sensitivity analysis and uncertainty analysis. After the 50000 runs, the threshold of sum of squared errors (SSE) for the objective function was taken to 2, and then the following plots were obtained.


### 4.3.1. Dotty Plots of Objective Function

The dotty plots shown in Figure1 is the result of rejection threshold of sum of squared errors (SSE) set to 2. The “surfaces” of the dotty plot for a, b, and c parameters have a clearly defined minimum, so they can be considered to be well identified. Followed by E
a, however, the A parameter is difficult to be defined.Figure 1
Dotty plots of objective function (SSE).


### 4.3.2. Posteriori Parameter Distribution

The initial simulation population was uniformly distributed, and the resulting distribution valuesD are plotted as bars as shown in Figure 2. In plots, the range of every parameter is divided into 20 containers of equal width. The likelihoods in each container are added up and divided by the sum of parameter values within the container. The regions of higher bars in Figure 2 show that their corresponding parameter ranges have better model performance. So, we can obtain the better region of parameters aimed at the model from Figure 2.Figure 2
Posteriori parameter distribution.


### 4.3.3. Regional Sensitivity Analysis Plot

The parameter sets are then split into 10 groups. For each group the likelihoods are normalized by dividing by their total, and the cumulative frequency distribution is calculated and plotted. If the model performance is sensitive to a particular parameter, there will be a large difference between the cumulative frequency distributions. Figure3 shows that parameter c is most sensitive, followed by parameters a, b, and E
a. And parameter A is insensitive.Figure 3
Regional sensitivity analysis plot.


### 4.3.4. 3D Surface Plot of Relationship between the Two Parameters

In order to examine parameter interaction, a Monte Carlo simulation was run with just 2 parameters varying and the other parameters fixed for MCAT. The plots in Figure4 show the parameter surfaces of every two parameters. As seen in Figure 4, we can note that there is a relatively well-defined region of low SSE, which indicates better combination state of the parameters. Therefore, it can provide the basis to judge the results of parameter identification.Figure 4
3D surface plot of every two parameters.
(a)
A
-
E
a
(b)
A-a
(c)
A-b
(d)
A-c
(e)
E
a-a
(f)
E
a-b
(g)
E
a-c
(h)
a-b
(i)
a-c
(j)
b-c


### 4.3.5. GLUE Variable Uncertainty

Figure5 shows the cumulative probability distribution and probability density function for the output variable calculated using a selected objective (transformed to likelihood). It also shows the uncertainty in peak output for the above example. The result shows that there is a probability of 0.05 (5%) that the peak output will be above 0.999 or below 0.865.Figure 5
GLUE variable uncertainty plot.The above results can help us to identify the important model input parameters, identify the main factors of uncertainty for model output, improve the model structure, and judge the rationality of identification results.


## 4.3.1. Dotty Plots of Objective Function

The dotty plots shown in Figure1 is the result of rejection threshold of sum of squared errors (SSE) set to 2. The “surfaces” of the dotty plot for a, b, and c parameters have a clearly defined minimum, so they can be considered to be well identified. Followed by E
a, however, the A parameter is difficult to be defined.Figure 1
Dotty plots of objective function (SSE).


## 4.3.2. Posteriori Parameter Distribution

The initial simulation population was uniformly distributed, and the resulting distribution valuesD are plotted as bars as shown in Figure 2. In plots, the range of every parameter is divided into 20 containers of equal width. The likelihoods in each container are added up and divided by the sum of parameter values within the container. The regions of higher bars in Figure 2 show that their corresponding parameter ranges have better model performance. So, we can obtain the better region of parameters aimed at the model from Figure 2.Figure 2
Posteriori parameter distribution.


## 4.3.3. Regional Sensitivity Analysis Plot

The parameter sets are then split into 10 groups. For each group the likelihoods are normalized by dividing by their total, and the cumulative frequency distribution is calculated and plotted. If the model performance is sensitive to a particular parameter, there will be a large difference between the cumulative frequency distributions. Figure3 shows that parameter c is most sensitive, followed by parameters a, b, and E
a. And parameter A is insensitive.Figure 3
Regional sensitivity analysis plot.


## 4.3.4. 3D Surface Plot of Relationship between the Two Parameters

In order to examine parameter interaction, a Monte Carlo simulation was run with just 2 parameters varying and the other parameters fixed for MCAT. The plots in Figure4 show the parameter surfaces of every two parameters. As seen in Figure 4, we can note that there is a relatively well-defined region of low SSE, which indicates better combination state of the parameters. Therefore, it can provide the basis to judge the results of parameter identification.Figure 4
3D surface plot of every two parameters.
(a)
A
-
E
a
(b)
A-a
(c)
A-b
(d)
A-c
(e)
E
a-a
(f)
E
a-b
(g)
E
a-c
(h)
a-b
(i)
a-c
(j)
b-c


## 4.3.5. GLUE Variable Uncertainty

Figure5 shows the cumulative probability distribution and probability density function for the output variable calculated using a selected objective (transformed to likelihood). It also shows the uncertainty in peak output for the above example. The result shows that there is a probability of 0.05 (5%) that the peak output will be above 0.999 or below 0.865.Figure 5
GLUE variable uncertainty plot.The above results can help us to identify the important model input parameters, identify the main factors of uncertainty for model output, improve the model structure, and judge the rationality of identification results.


## 4.4. Parameter Estimation

62 groups of experimental data including temperatureT, pressure P, residence time τ, [
2
CP
], [
O
2
], [
H
2
O
], and 2-CP conversion X provided in the literature [28] are used as test samples to verify the efficiency of the proposed algorithms.In the experiments, the population of DS is set to 50, the total iterative number is 500, scale factor is produced by using the gamma distribution, others are set with default parameters, and the parameter ranges of are[
0,200
], [
30000,50000
], [
0,2
], [
0,2
], and [
0,2
], respectively. Taking into account the situation in which DS has shorter running time, we increased an experiment where the total iterative number is 10000.As a comparison, the population of QPSO is set to 50; the total iterative number is 500;α decreases linearly from 1 to 0.5; the population of bees in ABC is set to 50; the total iterative number is 500; the limit value is 100; the ranges of parameters are the same as the previous ones.In order to prevent influence of random factor, the above algorithms are executed in 20 independent runs and the best, average, and the worst values of the results are recorded. The optimization results of the above algorithms are listed in Table1 in detail. As a comparison, the best results from some other algorithms are as follows: 0.2225 (chaos GA), 0.2494 (nonlinear regression), 0.2177 (adaptive differential evolution), 0.2177 (hybrid ant colony system), 0.2753 (RNA-GA, [13]), and 0.2551 (MDNA-GA, [13]).Table 1
Objective function values obtained by different algorithms.
Algorithms
Best
Average
Worst
DE [10]
0.217685
0.217725
0.218037
PSO [10]
0.217685
0.217819
0.218733
Hierarchical differential evolution [10]
0.217685
0.217685
0.217685
ABC [12]
0.233919
0.292802
0.421080
Hybrid ABC (HABC) [12]
0.217685
0.217685
0.217685
QPSO [29]
0.2176855
0.2177133
0.2178465
DS(1)
0.222266
0.252299
0.287777
DS(2)
0.2176850
0.2176855
0.217689
DS(lhs)
0.222360
0.248886
0.308399
Improved DS(1)
0.2176850
0.2176850
0.2176850
Improved DS(2)
0.2176850
0.2176850
0.2176850
Note: DS(1) represents the iteration number which is 500; DS(2) represents the iteration number which is 10000; DS(lhs) represents the initial population which is generated by Latin hypercube sampling; Improved DS(1) represents the simplex method which is used according to the sequential order; Improved DS(2) represents the other hybrid method, whose total iteration number is 5.Table1 shows that the best objective function values achieved by the algorithms (DE, PSO, HDE, HABC, QPSO, DS(2), Improved DS(1), and Improved DS(2)) are very similar, while the average and worst values obtained by HDE, hybrid ABC, DS(2), and the two improved DS algorithms outperform those of DE, PSO, and ABC. Furthermore, the basic DS algorithm can be used successfully with a large number of runs, and DS(lhs) does not show better optimization results. The obtained results also show that the two improvement algorithms compared with simplex obtain better performance, which has the advantages of strong global search ability and local search ability, can complete the optimization task outstandingly. The optimum solution (A
,
E
a
,
a
,
b
,
c) is 63.5468, 45626, 0.8081, 0.4444, and 0.3239, respectively.


## 5. Conclusions

In this paper, a new population-based heuristic optimization algorithm (DS) and its improved algorithms are applied to estimate the kinetic parameter of a chemical model. Main conclusions are as follows.(1)
Sensitivity analysis using MCAT showed that parameterc is most sensitive, followed by parameters a, b, and E
a, and A is insensitive. On the degree of difficulty of identification, a, b, and c are well identified, E
a is the second; however, A is difficult to be defined. Furthermore, by comparing, the optimization results of the intelligent algorithms are consistent with the distribution feature from parameter distribution plots and are corresponding to those regions which are indicated in the relationship plots of parameters. So, the results from MCAT can effectively help us to identify the important model input parameters, identify the main factors of uncertainty for model output, and judge the rationality of identification results.
(2)
DS algorithm has been compared with some widely used different optimization algorithms (DE, PSO, HDE, ABC, HABC, and QPSO) in solving the same problem. The results show that, with the large number of runs, the basic DS algorithm can be used successfully for solution of the above problem. But then, aimed at the shortcoming of the basic DS, the paper proposed two improvement measures. The obtained results show that the first improvement algorithm (DS-lhs) cannot effectively obtain better performance, and the second improvement algorithm (hybrid method), which has the advantages of strong global search ability and local search ability, can complete the optimization task outstandingly.
(3)
The inverse problem has nonunique characteristic, also known as the ill-posed problem of inversion theory. The problem can be solved by sensitivity analysis and uncertainty analysis. Monte Carlo method is often used to make sensitivity analysis and uncertainty analysis. Recently, the emerging swarm intelligence algorithms, which are based on the population evolution, are similar to the Monte Carlo method in the running mechanism. Therefore, combining with uncertainty analysis technology, the swarm intelligence optimization algorithms have wide applications in parameter identification.



---

*Source: 313105-2015-08-11.xml*