Steady Modeling for an Ammonia Synthesis Reactor Based on a Novel CDEAS-LS-SVM Model

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 thenet value ofammonia. Themodeling method has fastconvergence 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.


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 2 Mathematical Problems in Engineering 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.

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 Q1 and Q2 to go to the converter.The first cold shot Q1 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 Q1 from the bottom of reactor goes through the preheater and is heated by the counter-current flowing reacted gas from waste heat boiler.Q1 gas is divided into 4 cold quench gas (q1, q2, q3, and q4) and Q2 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 N 2 , H 2 , NH 3 , 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.

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  min and  max :  0 , =  0 min, + rand (0, 1) ( 0 max, −  0 min, )  = 1, 2, . . ., . (2) (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 1  , 2  , 3  , 4  , and 5  are mutually exclusive integers randomly generated within the range [1, NP] which should not be . is the mutation factor for scaling the difference vector, usually bounded in [0, 2].  best 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. rand is an integer chosen randomly in [1, 𝐷].To ensure that the trial vector (   ) differs from its corresponding individual (   ) by at least one dimension,  =  rand 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.

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][17][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.

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 [0, 1],  ∈ {0.1×},  = 1, 2, . . ., 10 and   ∈ {0.1×},  = 1, 2, . . ., 10.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 0 ≤   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; best  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 rand 1 and rand 2 are two random values which are uniformly distributed in [0, 1].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 {0.2, 0.4, 0.6, 0.8, 1.0}, 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 0 ≤  2 < 1 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; best  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 rand 3 and rand 4 are two random values which are uniformly distributed in [0, 1].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.

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.[19].In this paper, 30% of the individuals in the belief space are replaced by the good ones in population space.

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.

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.

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

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  = 0.3 and   = 0.5; for CDEAS, the size of both two subpopulations is 25; the initial  and   are randomly selected in (0, 1) and the initial mutation strategy is DE/rand/1; the interval information exchanges between the two subpopulations  is 50 generations; the thresholds   =    = 0.5 and  1 =  2 = 0.1.
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  2 ,  4 ,  7 ,  9 ,  11 ,  13 , and  15 is shifted.Functions  1 ∼ 7 are unimodal and functions  8 ∼ 18 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  2 and  7 in all trials, and the success rate reached 100% of functions  1 ,  2 ,  3 ,  4 ,  6 ,  7 , and  18 .For most of the test functions, the success The convergence figures of CDEAS comparing with original DE for 18 instances are listed as Figure 6.
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.

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.

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 [ × 1] vector of ones,  is the transpose of a matrix or vector,  is a weight vector,  0 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.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.

Results and Discussion
Operational parameters such as  H 2 ,  CH 4 , and   were collected and acquired from plant DCS from the year 2011-2012.In addition, data on the inlet ammonia concentration of recycle gas  NH 3 were simulated by mechanism and soft sensor model [20].
The extreme values are eliminated from the data using the 3 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 crossvalidation.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.

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.

Figure 4 :
Figure 4: Relationship between pheromone and ant paths of ,   .

Figure 5 :
Figure 5: Relationship between and ant paths of mutation strategy.

Table 1 :
The comparison results of the CDEAS algorithm and original DE algorithm.

Table 1 :
Continued.CDEAS gets very close to the global optimum in some other functions  1 ,  3 ,  4 ,  6 , and  18 .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  1 ,  3 ,  4 ,  5 ,  6 ,  12 ,  13 ,  14 ,  15 , and  18 although the computing time of CDEAS is longer than that of original DE because of its complexity.

Table 2 :
Auxiliary variables of model of net value of Ammonia.
∘ C

Table 3 :
The comparisons of training error and testing error of LS-SVM.