Abstract
A steady-state mathematical model is built in order to represent plant behavior under stationary operating conditions. A novel modeling using LS-SVR based on Cultural Differential Evolution with Ant Search is proposed. LS-SVM is adopted to establish the model of the net value of ammonia. The modeling method has fast convergence speed and good global adaptability for identification of the ammonia synthesis process. The LS-SVR model was established using the above-mentioned method. Simulation results verify the validity of the method.
1. Introduction
Ammonia is one of the important chemicals that has innumerable uses in a wide range of areas, that is, explosive materials, pharmaceuticals, polymers, acids and coolers, particularly in synthetic fertilizers. It is produced worldwide on a large scale with capacities extending to about 159 million tons at 2010. Generally, the average energy consumption of ammonia production per ton is 1900 KG of standard coal in China, which is much higher than the advanced standard of 1570 KG around the world. At the same time, the haze and particulate matter 2.5 has been serious exceeded in big cities in China at recent years, and one of the important reasons is the emission of coal chemical factories. Thus, an economic potential exists in energy consumption of the ammonia synthesis as prices of energy rise and reduce the ammonia synthesis pollution to protect the environment. Ammonia synthesis process has the characteristics of nonlinearity, strong coupling, large time-delay and great inertia load, and so forth. Steady-state operation-optimization can be a reliable technique for output improvement and energy reduction without changing any devices.
The optimization of ammonia synthesis process highly relies on the accurate system model. To establish an appropriate mathematical model of ammonia synthesis process is a principal problem of operation optimization. It has received considerable attention since last century. Heterogeneous simulation models imitating different types of ammonia synthesis reactors have been developed for design, optimization and control [1]. Elnashaie et al. [2] studied the optimization of an ammonia synthesis reactor which has three adiabatic beds. The optimal temperature profile was obtained using the orthogonal collocation method in the paper. Pedemera et al. [3] studied the steady state analysis and optimization of a radial-flow ammonia synthesis reactor.
The above study indicated that both the productive capacity and the stability of the ammonia reactor are influenced by the cold quench and the feed temperature significantly. Babu and Angira [4] described the simulation and optimization design of an auto-thermal ammonia synthesis reactor using Quasi-Newton and NAG subroutine method. The optimal temperature trajectory along the reactor and optimal flows throughput 3.3% additional ammonia production. Sadeghi and Kavianiboroujeni [1] evaluated the process behavior of an industrial ammonia synthesis reactorby one-dimensional model and two-dimensional model; genetic algorithm (GA) was applied to optimize the reactor performance in varying its quench flows. From The above literatures we can find that most models are built based on thermodynamic, kinetic and mass equilibria calculations. It is very difficult to simulate the specific internal mechanism because a lot of parameters are unknown in real industrial process.
In order to achieve the required accuracy of the model, some researches focus on the novel modeling methods combining some heuristic methods such as ANN (Artificial Neural Network), LS-SVM (Least Squares Support Vector Machine) with Evolutionary Algorithm, for example, genetic algorithm, ant colony optimization (ACO), particle swarm optimization (PSO), differential evolution (DE), and so forth. DE is one of the most popular algorithms for this problem and has been applied in many fields. Sacco and Hendersonb [5] introduced a variant of the differential evolution algorithm with a new mutation operator based on a topographical heuristic, and used it to solve the nuclear reactor core design optimization problem. Rout et al. [6] proposed a simple but promising hybrid prediction model by suitably combining an adaptive autoregressive moving average architecture and differential evolution for forecasting of exchange rates. Ozcan et al. [7] carried out the cost optimization of an air cooling system by using Lagrange multipliers method, differential evolution algorithm and particle swarm optimization for various temperatures and mass flow rates. The results showed that the method gives high accuracy results within a short time interval. Zhang et al. [8] proposed a hybrid differential evolution algorithm for the job shop scheduling problem with random processing times under the objective of minimizing the expected total tardiness. Arya and Choube [9] described a methodology for allocating repair time and failure rates to segments of a meshed distribution system using differential evolution technique. Xu et al. [10] proposed a model of ammonia conversion rate by LS-SVM and a hybrid algorithm of PSO and DE is described to identify the hyper-parameters of LS-SVM.
To describe the relationship between net value of ammonia in ammonia synthesis reactor and the key operational parameters, least squares support vector machine is employed to build the structure of the relationship model, in which a novel algorithm called CDEAS is proposed to identify the parameters. The experiment results showed that the proposed CDEAS-LS-SVM optimizing model is very effective of being used to obtain the optimal operational parameters of ammonia synthesis converter.
The remaining of the paper is organized as follows. Section 2 describes the ammonia synthesis production process. Section 3 proposes a novel Cultural Differential Evolution with Ant Colony Search (CDEAS) algorithm. Section 4 constructs a model using LS-SVM based on the proposed CDEAS algorithm. Section 5 presents the experiments and computational results and discussion. Finally, Section 6 summarizes the above results and presents several problems which remain to be solved.
2. Ammonia Synthesis Production Process
A normal ammonia production flow chart includes the synthesis gas production, purification, gas compression, and ammonia synthesis. Ammonia synthesis loop is one of the most critical units in the entire process. The system has been realized by LuHua Inc., a medium fertilizers factory of YanKuang Group, China.
Figure 1 represents a flow sheet for the ammonia synthesis process. The ammonia synthesis reactor is a one-axial flow and two-radial flow three-bed quench-type unit [11]. Hydrogen-nitrogen mixture is reacted in the catalyst bed under high temperature and pressure. The temperature in the reactor is sustained by the heat of reaction because the reaction is exothermic [1]. The reaction of ammonia synthesis process contains
The reaction is limited by the unfavorable position of the chemical equilibrium and by the low activity of the promoted iron catalysts with high pressure and temperature [12]. In general, no more than 20% of the synthesis gas is converted into ammonia per pass even at high pressure of 30 MPa [12]. As the ammonia reaction is exothermic, it is necessary for removing the heat generated in the catalyst bed by the progress of the reaction to obtain a reasonable overall conversion rate as same as to protect the life of the catalyst [13]. The mixture gas from the condenser is divided into two parts and to go to the converter. The first cold shot is recirculated to the annular space between the outer shell reactor and catalyst bed from the top to the bottom to refrigerate the shell and remove the heat released by the reaction. Then the gas from the bottom of reactor goes through the preheater and is heated by the counter-current flowing reacted gas from waste heat boiler. gas is divided into 4 cold quench gas (, , , and ) and gas for mixing with the gas between consecutive catalyst beds to quench the hot spots before entry to the subsequent catalyst beds. The hot spot temperatures (TIRA705, TIRA712N, and TIRA714) represent the highest reaction temperatures at each stage of the catalyst bed.
Figure 2 represents the ammonia synthesis unit. The reacted gas including N2, H2, NH3, and inert gas after reactor passes through the waste heat boiler. Then it goes through the preheater and the water cooler to be further cooled. Part of the ammonia is condensed and separated by ammonia separator I. Inert gas from the ammonia synthesis loop are ejected by purge gas from separator to prevent accumulation of inert gas in the system. The fresh feed gas is produced by the Texaco coal gasification air separation section, a process that converts the Coal Water Slurry into synthesis gas for ammonia. The fresh gas consists of hydrogen and nitrogen in stoichiometric proportions of 3 : 1 approximately and mixes with small amounts of argon and methane. The fresh gas which passes compressor is compounded with the recycle gas which comes from the circulator, and then the mixture goes through oil separator and condenser. Mixture gas is further cooled by liquid ammonia and goes through ammonia separator II to separate the partial liquid ammonia, and then it goes out with very few ammonia. The liquid ammonia from ammonia separator I and separator II flows to the liquid ammonia jar. Mixture is heated in ammonia condenser to about 25°C and flows to the reactor and the whole cycle starts again.
3. Proposed Cultural Differential Evolution with Ant Search Algorithm
3.1. Differential Evolution Algorithm
Evolutionary Algorithms, which are inspired by the evolution of species, have been adopted to solve a wide range of optimization problems successfully in different fields. The primary advantage of Evolutionary Algorithms is that they just require the objective function values, while properties such as differentiability and continuity are not necessary [14].
Differential evolution, proposed by Storn and Price, is a fast and simple population based stochastic search technique [15]. DE employs mutation, crossover, and selection operations. It focuses on differential vectors of individuals with the characteristics of simple structure and rapid convergence. The detailed procedure of DE is presented below.
(1) Initialization. In a -dimension space, NP parameter vectors so-called individuals cover the entire search space by uniformly randomizing the initial individuals within the search space constrained by the minimum and maximum parameter bounds and :
(2) Mutation. DE employs the mutation operation to produce a mutant vector called target vector corresponding to each individual after initialization. In iteration , the mutant vector of individual can be generated according to certain mutation strategies. Equations (3)–(7) indicate the most frequent mutation strategies version, respectively: where , , , , and are mutually exclusive integers randomly generated within the range which should not be . is the mutation factor for scaling the difference vector, usually bounded in . is the best individual with the best fitness value at generation in the population.
(3) Crossover. The individual and mutant vector are hybridized to compose the trial vector after mutation operation. The binomial crossover is adopted by the DE in the paper, which is defined as where rand is a random number between in 0 and 1 distributed uniformly. The crossover factor is a probability rate within the range 0 and 1, which influences the tradeoff between the ability of exploration and exploitation. is an integer chosen randomly in . To ensure that the trial vector () differs from its corresponding individual () by at least one dimension, is recommended.
(4) Selection. When a newly generated trial vector exceeds its corresponding upper and lower bounds, it is reinitialized within the presetting range uniformly and randomly. Then the trial individual is compared with the individual , and the one with better fitness is selected as the new individual in the next iteration:
(5) Termination. All above three evolutionary operations continue until termination criterion is achieved, such as the evolution reaching the maximum/minimum of function evaluations.
As an effective and powerful random optimization method, DE has been successfully used to solve real world problems in diverse fields both unconstrained and constrained optimization problems.
3.2. Cultural Differential Evolution with Ant Search
As we mentioned in Section 3.1, mutation factor , mutation strategies, and crossover factor have great influence on the balance of DE’s exploration and exploitation ability. decides the amplification of differential variation; is used to control the possibility of the crossover operation; mutation strategies have great influence on the results of mutation operation. In some literatures , , and mutation strategies are defined in advance or varied by some specific regulations. But the factors , , and strategies are very difficult to choose since the prior knowledge is absent. Therefore, Ant Colony Search is used to search the suitable combination of , , and mutation strategies adaptively to accelerate the global search. Some researchers have found an inevitable relationship between the parameters (, , and mutation strategies) and the optimization results of DE [16–18]. However, the approaches above are not applying the most suitable , , and mutation strategies simultaneously.
In this paper, based on the theory of Cultural Algorithm and Ant Colony Optimization (ACO), an improved Cultural Differential Algorithm incorporation with Ant Colony Search is presented. In order to accelerate searching out the global solution, the Ant Colony Search is used to search the optimal combination of and in subpopulation 1 as well as mutation strategy in subpopulation 2. The framework of Cultural Differential Evolution with Ant Search is briefly described in Figure 3.
3.2.1. Population Space
The population space is divided into two parts: subpopulation 1 and subpopulation 2. The two subpopulations contain equal number of the individuals.
In subpopulation 1, the individual is set as ant at each generation. and are defined to be the values between , , and , . Each of the ants chooses a combination of and according to the information which is calculated by the fitness function of ants. During search process, the information gathered by the ants is preserved in the pheromone trails . By exchanging information according to pheromone, the ants cooperate with each other to choose appropriate combination of and . Then ant colony renews the pheromone trails of all ants.
Then, the pheromone trail is updated in the following equation: where means the pheromone trail evaporation rate, , ; 1st parameter represents and 2nd parameter represents ; is the quantity of the pheromone trail of ant , where is the ant group that chooses th value as the selection of th parameter; denotes the best individual of ant colony till th generation.
In order to prevent the ants from being limited to one ant path and improve the possibility of choosing other paths considerably, the probability of each ant chooses th value of th parameter ( and ) in Figure 4 is set by
Figure 4 illustrates the relationship between pheromone matrix and ant path of and , where is a constant which is defined as selection parameter and and are two random values which are uniformly distributed in . Selection of the values of and depends on the pheromone of each path. According to the performance of all the individuals, the individual is chosen by the most appropriate combination of and in each generation.
In subpopulation 2, the individual is set as ant at each generation. Mutation strategies which are listed at (3)–(7) are defined to be of the values , respectively. For example, 0.2 means the first mutation strategy equation (3) is selected. Each of the ants chooses a mutation strategy according to the information which is calculated by the fitness function of ants. During search process, the information gathered by the ants is preserved in the pheromone trails . By exchanging information according to pheromone, the ants cooperate with each other to choose appropriate mutation strategy. Then ant colony renews the pheromone trails of all ants.
Then, the pheromone trail is updated in the following equation: where means the pheromone trail evaporation rate and is the quantity of the pheromone trail of ant , where is the ant group that chooses th value as the selection of parameter; denotes the best individual of ant colony till th generation.
In order to prevent the ants from being limited to one ant path and improve the possibility of choosing other paths considerably, the probability of each ant choosing th value of th parameter (mutation strategies) is set by where is a constant which is defined as selection parameter and and are two random values which are uniformly distributed in . Selection of the values of mutation strategies depends on the pheromone of each path. According to the performance of all the individuals, the individual is chosen by the most appropriate combination of mutation strategies in each generation.
Figure 5 illustrates the relationship between pheromone matrix and ant path of mutation strategies.
3.2.2. Belief Space
In our approach, the belief space is divided into two knowledge sources, situational knowledge and normative knowledge.
Situational knowledge consists of the global best exemplar which is found along the searching process and provides guidance for individuals of population space. The update of the situational knowledge is done if the best individual found in the current populations space is better than .
The normative knowledge contains the intervals that decide the individuals of population space where to move. and are the lower and upper bounds of the search range in population space. and are the value of the fitness function associated with that bound. If the and are updated, the and must be updated too.
and are set by
3.2.3. Acceptance Function
Acceptance function controls the amount of good individuals which impact on the update of belief space [19]. In this paper, 30% of the individuals in the belief space are replaced by the good ones in population space.
3.2.4. Influence Function
In the CDEAS, situational knowledge and normative knowledge are involved to influence each individual in the population space, and then population space is updated.
The individuals in population space are updated in the following equation: where is a constant of 0.2.
3.2.5. Knowledge Exchange
After steps, the and of subpopulation 2 are replaced by the suitable and calculated by subpopulation 1 and the mutation strategy of subpopulation 1 is displaced by the suitable mutation strategy calculated by subpopulation 2 simultaneously. So the and and mutation strategy are varying in the two subpopulations to enable the individuals to converge globally and fast.
3.2.6. Procedure of CDEAS
The procedure of CDEAS is proposed as follows.
Step 1. Initialize the population spaces and the belief spaces; the population space is divided into subpopulation 1 and subpopulation 2.
Step 2. Evaluate each individual’s fitness.
Step 3. To find the proper , , and mutation strategy, the Ant Colony Search strategy is used in subpopulation 1 and subpopulation 2, respectively.
Step 4. According to acceptance function, choose good individuals from subpopulation 1 and subpopulation 2, and then update the normative knowledge and situational knowledge.
Step 5. Adopt the normative knowledge and situational knowledge to influence each individual in population space through the influence functions, and generate two corresponding subpopulations.
Step 6. Select individuals from subpopulation 1 and subpopulation 2, and update the belief spaces including the two knowledge sources for the next generation.
Step 7. If the algorithm reaches the given times, exchange the knowledge of , , and mutation strategy between subpopulation 1 and subpopulation 2; otherwise, go to Step 8.
Step 8. If the stop criteria are achieved, terminate the iteration; otherwise, go back to Step 2.
3.3. Simulation Results of CDEAS
The proposed CDEAS algorithm is compared with original DE algorithm. To get the average performance of the CDEAS algorithm 30 runs on each problem instance were performed and the solution quality was averaged. The parameters of CDEAS and original DE algorithm are set as follows: the maximum evolution generation is 2000; the size of the population is 50; for original DE algorithm and ; for CDEAS, the size of both two subpopulations is 25; the initial and are randomly selected in and the initial mutation strategy is DE/rand/1; the interval information exchanges between the two subpopulations is 50 generations; the thresholds and .
To illustrate the effectiveness and performance of CDEAS algorithm for optimization problems, a set of 18 representative benchmark functions which were listed in the appendix were employed to evaluate them in comparison with original DE. The test problems are heterogeneous, nonlinear, and numerical benchmark functions and the global optimum for , , , , , , and is shifted. Functions ~ are unimodal and functions ~ are multimodal. The detailed principle of functions is presented in [11]. The comparisons results of CDEAS and original DE algorithm are shown in Table 4 of the appendix. The experimental results of original DE and CDEAS algorithm on each function are listed in Table 1. Mean, best, worst, std., success rate, time represent the mean minimum, best minimum, worst minimum, the standard deviation of minimum, the success rate, and the average computing time in 30 trials, respectively.
From simulation results of Table 1 we can obtain that CDEAS reached the global optimum of and in all trials, and the success rate reached 100% of functions , , , , , , and . For most of the test functions, the success rate of CDEAS is higher in comparison with original DE. Moreover, CDEAS gets very close to the global optimum in some other functions , , , , and . It also presents that the mean minimum, best minimum, worst minimum, the standard deviation of minimum, and the success rate of CDEAS algorithm are clearly better than the original DE for functions , , , , , , , , , and although the computing time of CDEAS is longer than that of original DE because of its complexity.
The convergence figures of CDEAS comparing with original DE for 18 instances are listed as Figure 6.
From Figure 6 one can observe that the convergence speed of CDEAS is faster than original DE for , , , , , , , , , , , and .
All these comparisons of CDEAS with original DE algorithm have shown that CDEAS is a competitive algorithm to solve all the unimodal function problems and most of the multimodal function optimization problems listed above. As shown in the descriptions and all the illustrations before, CDEAS is efficacious on those typical function optimizations.
4. Model of Net Value of Ammonia Using CDEAS-LS-SVM
4.1. Auxiliary Variables Selection of the Model
There are some process variables which have the greatest influence on the net value of ammina, such as system pressure, recycle gas flow rate, feed composition (H/N ratio), ammonia and inert gas cencetration in the gas of reactor inlet, hot spot temperatures, and so forth. The relations between the process variables are coupling and the operational variables interact with each other.
The inlet ammonia concentration is an important process variable which is beneficial to operation-optimization but the device of online catharometer is very expensive. According to the mechanism and soft sensor model, a IIO-BP model was built to get the more accurate value of the inlet ammonia concentration [20]
From the analysis discussed above, some important variables have significant effects on the net value of ammonia. By discussion with experienced engineers and taking into consideration a priori knowledge about the process, the system pressure, recycle gas flow rate, the H/N ratio, hot-spot temperatures in the catalyst bed, and ammonia and methane concentration in the recycle gas are identified as the key auxiliary variables to model net value of ammonia which is listed in Table 2.
4.2. Modeling the Net Value of Ammonia Using CDEAS-LS-SVM
LS-SVM is an alternate formulation of SVM, which is proposed by Suykens. The e-insensitive loss function is replaced by a squared loss function, which constructs the Lagrange function by solving the problem linear Karush_Kuhn_Tucker (KKT) where is a vector of ones, is the transpose of a matrix or vector, is a weight vector, means the model offset, and is regression vector.
is Mercer kernel matrix, which is defined as where is defined by kernel function.
There are several kinds of kernel functions, such as hyperbolic tangent, polynomial, and Gaussian radial basis function (RBF) which are commonly used. Literatures have proved that RBF kernel function has strong generalization, so in this study RBF kernel was used: where and indicated different training samples, is the kernel width parameter.
As we can see from (19)~(21), only two parameters are needed for LS-SVM. It makes LS-SVM problem computationally easier than SVR problem.
Grid search is a commonly used method to select the parameters of LS-SVM, but it is time-consuming and inefficient. CDEAS algorithm has strong search capabilities, and the algorithm is simple and easy to implement. Therefore, this paper proposes the CDEAS algorithm to calculate the best parameters of LS-SVM.
5. Results and Discussion
Operational parameters such as , , and were collected and acquired from plant DCS from the year 2011-2012. In addition, data on the inlet ammonia concentration of recycle gas were simulated by mechanism and soft sensor model [20].
The extreme values are eliminated from the data using the criterion. After the smoothing and normalization, each data group is divided into 2 parts: 223 groups of training samples which are used to train model while 90 groups of testing samples which are valuing the generalization of the model for identifying the parameters of the LS-SVM, the kernel width parameter, and the weight vector.
BP-NN, LS-SVM, and DE-LS-SVM are also used to model the net value of ammonia, respectively. BP-NN is a 13-15-1 three-layer network with back-propagation algorithm. LS-SVM gains the with grid-search and cross-validation. The parameter settings of CDEAS-LS-SVM are the same as those in the benchmark tests. Each model is run 30 times and the best value is shown in Table 3. Descriptive statistics of training results and testing results of model include the relative error, absolute error, and mean square error. The performance of the four models is compared as shown in Table 3. The training and testing results of four models are illustrated in Figure 7.
Despite the fact that the training error using BP-NN is smaller than that using CDEAS-LS-SVM, which is because BP-NN is overfitting to the training data, the mean square error (MSE) on training data using CDEAS-LS-SVM is reduced by 25.6% and 23.2% compared with LS-SVM and DE-LS-SVM, respectively. In comparison with the other models (BP-NN, LS-SVM, and DE-LS-SVM), testing error using CDEAS-LS-SVM model is reduced by 14.1% and 11.2%, respectively. The results indicate that the proposed CDEAS-LS-SVM model has a good tracking precision performance and guides production better.
6. Conclusion
In this paper, an optimizing model which describes the relationship between net value of ammonia and key operational parameters in ammonia synthesis has been proposed. Some representative benchmark functions were employed to evaluate the performance of a novel algorithm CDEAS. The obtained results show that CDEAS algorithm is efficacious for solving most of the optimization problems comparisons with original DE. Least squares support vector machine is used to build the model while CDEAS algorithm is employed to identify the parameters of LS-SVM. The simulation results indicated that CDEAS-LS-SVM is superior to other models (BP-NN, LS-SVM, and DE-LS-SVM) and meets the requirements of ammonia synthesis process. The CDEAS-LS-SVM optimizing model makes it a promising candidate for obtaining the optimal operational parameters of ammonia synthesis process and meets the maximum benefit of ammonia synthesis production.
Appendix
(1) Sphere function
(2) Shifted sphere function
(3) Schwefel’s Problem 1.2
(4) Shifted Schwefel’s Problem 1.2
(5) Rosenbrock’s function
(6) Schwefel’s Problem 1.2 with noise in fitness
(7) Shifted Schwefel’s Problem 1.2 with noise in fitness
(8) Ackley’s function
(9) Shifted Ackley’s function
(10) Griewank’s function
(11) Shifted Griewank’s function
(12) Rastrigin’s function
(13) Shifted Rastrigin’s function
(14) Noncontiguous Rastrigin’s function
(15) Shifted noncontiguous Rastrigin’s function
(16) Schwefel’s function
(17) Schwefel’s Problem 2.21
(18) Schwefel’s Problem 2.22
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors are grateful to the anonymous reviewers for giving us helpful suggestions. This work is supported by National Natural Science Foundation of China (Grant nos. 61174040 and 61104178) and Fundamental Research Funds for the Central Universities, Shanghai Commission of Science and Technology (Grant no. 12JC1403400).