Nonlinear Model Predictive Control for pH Neutralization Process based on SOMA Algorithm

In this work, the pH neutralization process is described by a neural network Wiener (NNW) model. A nonlinear Model Predictive Control (NMPC) is established for the considered process. The main difficulty that can be encountered in NMPC is solving the optimization problem at each sampling time to determine an optimal solution in finite time. The aim of this paper is the use of global optimization method to solve the NMPC minimization problem. Therefore, we propose in this work, to use the Self Organizing Migrating Algorithm (SOMA) to solve the presented optimization problem. This algorithm proves its efficiency to determine the optimal control sequence with a lower computation time. Then the NMPC is compared to adaptive PID controller, where we propose to use the SOMA algorithm to formulate the PID in order to determine the optimal parameters of the PID. The performances of the two controllers based on the SOMA algorithm are tested on the pH neutralization process. Keywords—Nonlinear model predictive control; optimization; SOMA algorithm; adaptive PID; pH neutralization process


I. INTRODUCTION
The pH neutralization process is characterized by a nonlinear behavior.The high nonlinear characteristic of this process makes the control of the pH a hard task.Since the necessity of maintaining the pH in a specific range value, many control strategies have been proposed.For this purpose, [1]- [3] designate a PID controller to control the pH value.An adaptive controller is developed in [4].[5] proposed a model algorithmc control strategy.Other works [6]- [8] developed the Model predictive control (MPC) which is a more advanced control strategy.The MPC control strategy presents the major advantage to efficiently handle nonlinearity and constraints imposed on system input and output [10].Indeed, for real processes, the control values are usually delimited by an upper and lower bound that should be respected.
The MPC is essentially based the choice of a suitable model and adequate optimization method to solve the minimization problem.Many structures was used to describe the nonlinear behavior of the pH neutralization process: NARX model [11],Fuzzy neural network model [12], Neural networks [13], Wiener model [8], [14]- [17].
Respecting the nonlinear nature of the model, the resulting NMPC minimization problem is nonlinear and nonconvex.The nonconvexity will complicate the implementation of the NMPC.Added to that, solving a nonlinear optimization problem is a hard task in terms of computation time and burden.To overcome these difficulties a variety of solutions were proposed in literature to avoid solving a nonlinear optimization problem and ensure global convergence at each sampling time.In [18], [19], the minimization problem is converted into a linear one in the case of a simple polynomial description of the nonlinear block.In this case, the nonlinearity is removed by considering the inverse of the polynomial function.As a result, the input control will be linearly expressed in the predicted output leading to a quadratic optimization problem.Also, the nonlinear optimization problem is formulated in [8], [15], [16] as a linear problem by performing a linearization of the nonlinear model.To avoid solving the nonlinear optimization problem and reduce the implementation complexity of the NMPC, [20], proposed to describe the nonlinear process by a set of uncertain linear models instead of one nonlinear model.We can note that all these works avoid to solve the nonconvex optimization problem regarding the difficulty of implementation and the high computation burden necessary at each sampling time.
To obtain good control accuracy the NMPC optimization problem must be solved as nonlinear [21].However, global optimization methods are generally not recommended since they are not able to ensure the real time feasibility of the NMPC.In fact, the sampling period of the process under consideration must be respected at each iteration when solving the NMPC problem.This constraint is difficult to be satisfied using global optimization method due to their slow convergence.Therefore, the main aim of this work is the use of an efficient algorithm to solve the NMPC optimization problem and ensure the realtime feasibility of the control algorithm.For this, we propose to use the Self Organizing Migrating algorithm (SOMA) in this work to solve the optimization problem.The SOMA is an evolutionary algorithm such as Genetic Algorithm (GA), Particle Swarm Optimization (PSO), Artificial Bee Colony algorithm (ABC) and so on.The SOMA algorithm was successfully applied to solve a variety of engineering problems.The most interesting are: control issues [22], antenna design [23], system identification [24], [25], Aircraft wing design and Synthesis of robot control program [26].This method presents the significant advantage of high convergence speed.
In This work, we will be interested to represent the pH neutralization process by a Wiener model.Wiener model belongs to block oriented models that are defined by a dynamic linear block followed by a nonlinear static one.Due to the specific structure of the pH neutralization process, the Wiener model appears able to well reproduce the behavior of the process.Many works used the Wiener model to describe the nonlinear dynamic of the process.
The different representations vary in the way that linear and nonlinear blocks are described.
In this paper, the NMPC is integrated with SOMA algorithm to solve the optimization problem.We prove, in this work, the ability of the SOMA algorithm to ensure good control performances with a low computation time despite the large prediction and control horizons.In the sequel, The NMPC strategy integrated with SOMA algorithm is compared to adaptive PID controller.We propose, in this work, to adjust PID parameters using SOMA algorithm.This paper is organized as follows: Section 2 describes the Wiener model.The identification of the considered process is detailed in Section 3. NMPC based on the Wiener model is presented in Section 4. The SOMA algorithm used to solve the minimization problem is described in Section 5. Simulation results are given in Section 6.

II. WIENER MODEL
Wiener model belongs to block oriented models.It is described by a linear dynamic block followed by a nonlinear static one as shown in Fig. 1.In this work, the linear dynamic block is described by a simple autoregressive model defined as: Where polynomials A and B are given by: A(q −1 ) = 1 + a 1 q −1 + a 2 q −2 +, ..., +a na q −na (2) n a and n b are respectively the orders of the two polynomials A and B. The nonlinear block is defined as: f (.) is a nonlinear function.
Different forms are used to describe the nonlinear block for Wiener model starting from the simple polynomial form [19], [28] to more complex description Neural network [19], [15], support vector machine [16].In this work the nonlinear static block of the Wiener model is given by a feed-forward neural network as adpoted in [15].

III. IDENTIFICATION OF THE WIENER MODEL
The aim of this paragraph is the determination of the parameters a i and b i of the linear block as well as the weights w i of the network that present the parameters of the nonlinear block.
The structure of the Wiener model is depicted in Fig. 2.

A. Identification of the Linear Block
Firstly, a small input signal is applied to the system to ensure linear perturbation of the nonlinear system [29].Then, the recursive least square algorithm is firstly used to identify the parameters of linear dynamic block.
The output of the linear block can be expressed as: ψ is the data vector defined as: θ is the parameter vector: The update of the parameter vector is ensured by the following equation: P is weighting matrix given by: The initial parameter vector is θ(0) = 0.
Once the parameters of the linear block are determined, the back-propagation algorithm is applied to train the feed-forward neural network such as in [29].

B. Identification of the Nonlinear Block
The structure of the nonlinear block is illustrated in Fig. 3.The output of the nonlinear block can be computed from Fig. 3 as: where j is the number of hidden nodes, f is a nonlinear activation function taken as the hyperbolic tangent function and e(k) is the output of the hidden layer defined as: The weights w i,j of the neural network are updated by minimizing the following criterion: where y(k) is the process output and y M (k) is the model output defined as: Applying the back-propagation training algorithm the optimization is carried out in order to minimize the criterion (12) with respect to the weights w i,j of the network.
The update equation of the different weights is defined as: where ∆w(k) is defined as: ∆w µ is the learning coefficient.
The Neural Network Wiener (NNW) model is used in this work to compute the j-step ahead predictions.Future prediction are used to define the cost function of the NMPC strategy.

IV. MODEL PREDICTIVE CONTROL DESIGN FOR WIENER MODEL
The basic idea of the MPC strategy is the use of the process model to compute the j-step ahead prediction over a prediction horizon N p .At each sampling time, a control sequence U = [u(k), u(k + 1), ..., u(k + N u − 1)] T is computed , where N u is the control horizon, by minimizing the cost function defined as the difference between the predicted output ŷ(k + i|k) and the future set-point y sp (k + i) defined by: Based on the receding control horizon, only the first element of the control sequence U is applied to the system.Then, the whole procedure will be repeated at the next sampling time and the prediction horizon will be shifted one step forward.
In order to get efficient control results, predictions must be accurately computed.It is so important to take into account the effects of system and model mismatch coming from modeling errors and unmeasured disturbance, a correction term as in the dynamical matrix control, is added to model output defined as [30]: where y M (k + i|k) is the model output and y(k) is the process output.The correction term is constant over the prediction horizon.Therefore, the model expression that will be used to compute predictions is defined as: Due to the nonlinear nature of the NNW model the resulting optimization problem will be nonlinear and nonconvex.The convergence of this optimization problem is not guaranteed and the algorithm may be trapped in a local minimum that will lead necessarily to suboptimal control performances.Solving a nonlinear optimization problem is a high time demanding task.Added to that, real process requires generally a large prediction and control horizons that will raise the computation time.Moreover, the sampling period of the process presents an additional constraint that should be respected when solving the NMPC optimization problem.Therefore, a fast convergence speed algorithm must be used to ensure the global convergence and the real-time feasibility of the control algorithm.Deterministic optimization methods are used in [9], [31] to solve the NMPC optimization problem.These methods are high time consuming also they can be only applied to systems requiring small prediction and control horizons.
Many research papers propose to use stochastic optimization methods to solve the minimization problem.The Artificial Bee Colony (ABC) algorithm is combined with the NMPC in [32] to solve the minimization problem.This algorithm proves its simplicity of implementation and reduced computational complexity.The genetic algorithm (GE) is used in [33], [34] to determine the optimal control sequence.The particle swarm optimization (PSO) algortihm is integrated in [35] with NMPC to solve the resulting optimization problem.In [36], the neural network is used to determine the solution of the minimization problem.In [37], the Nelder Mead algorithm was applied which leads to global solution by using far initialization.The simulations gave optimal results with least computation time for SISO and MIMO models.However, it remains a local optimization method.
We propose in this work to use the SOMA algorithm to solve the presented optimization problem.SOMA presents an effective, robust and simple global optimization method.

V. SELF ORGANIZING MIGRATING ALGORITHM
SOMA is a stochastic optimization method proposed first by Zelinka [25].SOMA is based on the social group of individual not on the philosophy of evolution.The classification of SOMA as evolutionary algorithm is explicated by the fact that the obtained results after a migration loop is the same as the result of one generation of evolutionary algorithm [26].The principle of SOMA algorithm can be summarized as a migration loops during which the position of each individual is enhanced in order to reach the leader position (individual with the best fitness).Each individual will be randomly initialized in the search space described by the upper and the lower bounds of the variables.At each migration loop, individuals are evaluated, the one that has the less fitness will be the leader, the rest individuals will cross a trajectory (pathlength) with step t in the direction of the leader.
Similar to other evolutionary algorithms, the operation of SOMA is ensured by some control parameters.The recommended ranges of these parameters are fixed based on great number of empirical tests and they are defined and given by [24], [26]: • Dim: It defines the problem dimension (number of decision variable).
• Population size (P opsize): It defines the number of individuals in population, recommended value ≥ 10.
• M igrations: It defines the maximum number of iterations (migration loops).It is the stopping criterion in SOMA recommended range [10, up to user].
• P athlength: It fixes how far the individual will stop its movement from the leader.If the pathlength=1 the individual will stop at the leader position, if pathlength=2 the individual will surpass the leader position by the same distance from the initial position, recommended value 3.
• step: It defines the step that uses individuals to cross the path.
• P RT : The P RT is used to determine the PRTvec.Individuals are allowed to change their position based on the P RT , recommended value 0.4.
SOMA is a population based algorithm.The initial population P is generated randomly in the search space defined by the lower x min and upper x max bound of the manipulated variables.So, P is defined as: P = {X 1 , X 2 , . . ., X P opsize }.The ith individual X i is defined by X i = {x i,1 , x i,2 , . . ., x i,Dim }, i = 1, . . ., P opsize, j = 1, . . ., Dim.
x i,j = x min,j + (x max,j − x min,j ).rand (20) Like other evolutionary algorithm, SOMA performs a set of stochastic evolutionary operators.These operators are defined by mutation and crossover.
• Mutation: This operator defines in evolutionary algorithm the diversity in the population.In SOMA mutation is applied differently.It is performed by the PRT vector noted PRTvec.The PRTvec aims to perturb the path of individuals randomly in order to ensure diversification among them [38].The perturbation in SOMA algorithm presents the mutation phase in the GA algorithm.Two values can be affected to this vector: 0 or 1 based on the SOMA PRT control parameter.Only individuals with PRTvec equal to 1 are allowed to change their positions.The PRTvec is created as follows: if rand < P RT P RT vec j = 1; else P RT vec j = 0; end • Crossover: Crossover operator means the creation of new individual during the search.Since in SOMA algorithm no new individual is generated, the crossover is defined in this case by generating a new best position of individual across the search space to reach the leader.At each migration loop individuals explore a set of positions when mapping the path, memorize the best found one and move to this position at the end of the path.At the next migration loop, individuals start from this position.The movement of individuals to reach the position of the leader is given by the following equation: where x k+1 i,j is the new individual position., x k i,j is the position of individual at iteration k, x k L,j is the leader position and t ∈ [0 : step : path length] The SOMA algorithm can be summarized as follows: 01 Choose the control parameters :PRT,step,pathlengthandmigration Generate the initial population P using equation (20), 02 Evaluate individuals of the population 03 Sort individuals and select the leader itetration = 1; while (iteration < migrations) f or i = 1 : P opsize Generate PRTvec using equation( 21) M ove each individual toward the leader using equation (22) Evaluate the individual at the new position new f itness < f itness(individual (i)) M ove individual(i) to the new best position.end if end f or sort individuals and select the new leader iteration = iteration + 1 end while

VI. SIMULATION RESULTS
The considered system is composed of the base stream q 1 , the acid stream q 2 and the buffer stream q 3 that are mixed in continuous stirred tank reactor.The dynamic model of the reactor is derived from the conservation equation and equilibrium relations [4].The dynamic of the pH process is given by two differential equations and a nonlinear one given by ( 23) and (24 ).
the different parameters of the reactor are listed in Table I.
In this section, we will consider the pH neutralization process.

A. Identification of the pH Neutralization Process
The pH neutralization process will be described by a NNW model, where the linear block is described by an autoregressive model and the nonlinear block is given by a multilayer feed-forward neural network with one hidden layer.The plant is excited using two different input signals to get the identification and the validation data as shown respectively in Fig. 4 and 5.The bounds on the input are 0 and 30.
A feed-forward neural network with one hidden layer and 3 nodes appears sufficient to describe the nonlinear block of the Wiener model [15].Fig. 6 and 7 illustrate the identification and the validation of the system.We can conclude that the neural Wiener model can reproduce the dynamic of the pH process with sufficient accuracy.This is well confirmed by the validation error depicted in Fig. 8 defined as the difference between the system and the model output y(k) and y M (k), respectively.

B. NMPC Control of the pH Neutralization Process
The NMPC strategy aim to maintain the pH value at a desired value.For this the neural Wiener model is used to compute the j-step-ahead predictions for the NMPC.The parameters of the NMPC are fixed as N p = 10, N u = 3, λ = 0.08, q min = 0 and q max = 30.
The recommended values of the parameters of the SOMA algorithm are fixed as: pathlength = 3, step = 0.31, P RT = 0.4, popsize = 10, the number of migration loops is 30.
We can note that the SOMA algorithm proves its efficiency to ensure good output tracking as depicted in Fig. 9.The performance offered by the SOMA algorithm outperforms those offered by the GBM method.This can be proved in terms of less overshoot and best control quality as illustrated in Fig. 9.In order to test the efficiency of the SOMA algorithm in presence of disturbance a constant v is added to the system output as: We can note from Fig. 10 the good ability of the SOMA algorithm to reject disturbance.This latter is rejected within 10 iterations that present a short time.
The NMPC strategy is compared to an adaptive PID controller.The Sum of square Error (SAE) defined by equation ( 25) is minimized in order to determine the proportional, integral and derivative coefficients of the PID controller: k p , k i and k d , respectively: Respecting the high nonlinear considered process, the SOMA algorithm is adopted to derive the coefficients of the adaptive PID controller at each iteration.
The control parameters of the SOMA algorithm are fixed for the PID controller as those of the NMPC.Simulation results for PID and NMPC are given in Fig. 11.The results show that both controllers ensure good output tracking in steady state for the first and second setpoint change.However, we can remark that the PID exhibits more overshoots as well as the deterioration of the control quality and consequently the output tracking in the last output change.This proves the superiority of the NMPC strategy to ensure good performances compared to PID.Computation time present an important performance index in control problem.Fig. 12 shows the computation time for both NMPC and PID.The present results confirm the good ability of NMPC to give better performances with the low computation time.In fact, the determination of optimal PID parameters on line at each sampling time using the SOMA algorithm requires more migration loops this will directly affect the computation time of the implementation of the PID controller.

VII. CONCLUSION
In this paper, a nonlinear model predictive control is combined with the SOMA algorithm to solve the NMPC optimization problem.The control performances of the SOMA algorithm are tested on a high nonlinear process.Control results prove the efficiency of this algorithm to ensure good output tracking and control accuracy with a low computation time.NMPC based on SOMA algorithm is compared to adaptive PID controller.We can conclude from simulation results that the NMPC outperforms the PID in terms of less overshoot and computation time.NMPC offers the best control results and the less computation time compared to PID.The difficulty of the determination of PID parameters for high nonlinear process limits the performance of this controller.

Fig. 9 .
Fig. 9. Output tracking using the SOMA algorithm and the GBM.

Fig. 10 .
Fig. 10.Output tracking in the presence of disturbance.

Fig. 12 .
Fig. 12. Computation time of the NMPC and PID