Constrained Nonlinear Model Predictive Control of an MMA Polymerization Process via Evolutionary Optimization

In this work, a nonlinear model predictive controller is developed for a batch polymerization process. The physical model of the process is parameterized along a desired trajectory resulting in a trajectory linearized piecewise model (a multiple linear model bank) and the parameters are identified for an experimental polymerization reactor. Then, a multiple model adaptive predictive controller is designed for thermal trajectory tracking of the MMA polymerization. The input control signal to the process is constrained by the maximum thermal power provided by the heaters. The constrained optimization in the model predictive controller is solved via genetic algorithms to minimize a DMC cost function in each sampling interval.


INTRODUCTION
Model predictive control (MPC), is a model based advanced process control (APC) technique that has been proved to be very successful in controlling highly complex dynamic systems. It naturally supports design for MIMO and time-delayed systems as well as state/input/output constraints. MPC is generally based on online optimization but in the case of unconstrained linear plants, closed form solutions can be derived analytically. However, when there are constraints over the control inputs (i.e. actuators) and/or process states, which is often the case, an online (i.e. real-time) constrained optimization problem has to be solved in each sampling interval, even if the plant model is linear and time invariant. This online optimization usually requires a high computational power; however, since chemical processes are typically of slow dynamics, such controllers have been designed and implemented on various chemical plants with great success. Moreover, due to recent advancements of computational hardware and software tools, the usage of MPC is rapidly expanding to other control domains including electrical machines, renewable energy, aerospace and automotive control systems.
In the past two decades, the effective control of polymerization processes control has been studied by many authors [1][2][3][4][5][6][7]. Polymerization kinetic is usually complex due to the nonlinearity of the process. Therefore, the control of the polymerization reactor has been staying a challenging task. Due to its great flexibility, a batch reactor is suitable to produce small amounts of special polymers and copolymers. The batch reactor is always dynamic by its nature. It is essential to have a suitable dynamic model of the process. Rafizadeh [1] presented a review on the proposed models and suggested an online estimation of some parameters. His model consists of the oil bath, electrical heaters, cooling water coil, and reactor. Peterson et al. [2] presented a non-linear predictive strategy for semi batch polymerization of MMA. Soroush and Kravaris [3] applied a Global Linearizing Control (GLC) method to control the reactor temperature. Performance of the GLC for tracking an optimum temperature trajectory was found to be suitable. DeSouza et al. [4] studied an expert neural network as an internal model in control of solution polymerization of vinyl estate. In their study, they compared their neural network control with a classic PID controller. Clarke-Pringle and MacGregor [5] studied the temperature control of a semi-batch industrial reactor. They suggested a coupled nonlinear strategy and extended Kalman filter method. Mutha et al. [6] suggested a non-linear model based control strategy, which includes a new estimator as well as Kalman filter. They conducted experiments in a small reactor for solution polymerization of MMA. Rho et al. [7] assumed a first order model plus dead time to pursue the control studies and estimated the parameters of this model by on line ARMAX model. Nonlinear predictive control of the batch reactor considered [8] and [9] via PCA and Wiener modeling approaches, respectively. When MPC is formulated as a state feedback controller, the full state information is required which must be provided using state estimators in nonlinear H2 (e.g. EKF) or nonlinear H∞ paradigms [10][11][12], [27][28]. Rafizadeh [13] designed a sequential linearization adaptive controller for the solution polymerization of methyl methacrylate in a Batch Reactor.
This paper presents a constrained model predictive control of an MMA reactor, based on the genetic algorithm optimization. A previously developed mechanistic model of the process was used. The model is a sequential piecewise linearization along a selected temperature trajectory. The piecewise linear model is used both for the plant output calculation through the prediction horizon and for the closed loop simulation of the controller, using a time-triggered switching mechanism. We are using an output feedback MPC, therefore, no state estimator is required, which is advantageous. The results of tracking the trajectory and eliminating noise and disturbances show a promising performance of the controller.

POLYMERIZATION MECHANISM
Methyl methacrylate normally is produced by a free radical, chain addition polymerization. Free radical polymerization consists of three main reactions: initiation, propagation and termination. Free radicals are formed by the decomposition of initiators. Once formed, these radicals propagate by reacting with surrounding monomers to produce long polymer chains; the active site being shifted to the end of the chain when a new monomer is added. During the propagation, millions of monomers are added to o P 1 radicals. During termination, due to reactions among free radicals, the concentration of radicals decreases. Termination is by combination or disproportionate reactions. With chain transfer reactions to monomer, initiator, solvent, or even polymer, the active free radicals are converted to dead polymer. Table 1 gives the basic free radical polymerization mechanism [14].
The free radical polymerization rate decreases due to reduction of monomer and initiator concentration. However, due to viscosity increase beyond a certain conversion there is a sudden increase in the polymerization rate. This effect is called Trommsdorff, gel, or auto-acceleration effect. For bulk polymerization of Methyl Methacrylate beyond the % 20 conversion, reaction rate and molecular weight suddenly increase. In high conversion, because of viscosity increase there is a reduction in termination reaction rate.

MATHEMATICAL MODELING OF POLYMERIZATION
The polymer production is accomplished by a reduction in volume of the mixture. The volumetric reduction factor is given by: The instantaneous volume of mixture is given by: The parameter β is defined as: During the free radical polymerization, the cage, glass, and gel effects occur. For the cage effect, the initiator efficiency factor is used. The CCS (Chiu, Carrat and Soong) model is used in this study to take into consideration the glass and the gel effects. Therefore, propagation rate constant, p k , is changing according to: Similarly, termination rate constant, t k , is given by θ and t θ are adjustable parameters related to propagation and termination rate constants, respectively. All other necessary parameters and constants for this model are given in [1], [14] and [15].
Long Chain Approximation (LCA) and Quasi Steady State Approximation (QSSA) are used in this study. Equations are highly nonlinear and, using Taylor expansion series, these equations were converted to linearized form. The linearized state space form is given by:  The molecular properties of the produced polymer are controlled by ensuring the reaction temperature is changing according to a desired reference trajectory. This is a tracking control problem which we are solving using MPC. Figure 1 shows the result of model validation [14]. As it is seen, the model is a good representative of the process. Equation (7) is converted to the transfer function for reaction temperature to input power:

Constrained Nonlinear Model Predictive Control of a Polymerization Process via Evolutionary Optimization
The result of sequential linearization is 131 transfer functions along the temperature profile. See [13] and [14] for a more detailed description of the MMA polymerization dynamic modeling.

Experimental Setup
A schematic representation of the experimental batch reactor setup is shown in Figure 2 [14]. The reactor is a Buchi type jacketed, cylindrical glass vessel. A multi-paddle agitator mixes the content. Two Resistance Temperature Detectors (RTDs) of PT100 type were used with accuracy of to measure the reactor temperature and the oil temperature in the oil bath. PT100 sensors provide good linearity in the measurement range and negligible drift. Methyl Methacrylate and Toluene were used as monomer and solvent, respectively. Benzoyl Peroxide (BPO) was used as the initiator. The heater, heats the oil circulating the oil bath, which is pumped into the reactor. Cool water is circulated in a coolant coil inside the oil bath through an electric on/off valve and acts as a safety feature to prevent the oil and consequently the reactor from overheating. The RTD outputs are converted into 0-10 VDC through a bridge and an instrument amplifier and are read by the data acquisition card A/D channels. The controller output is fed into a MOSFET-based power electronics switching circuit as PWM signals. The control command is applied to the MOSFETS through optical isolation, opto-couplers, which provide isolation from a three phase 220V, 50Hz power line. The maximum heater power available is 1000 W, which is a constraint on the control signal.

Model Predictive Control
Due to its high performance, model predictive control method has received a great deal of attention to control chemical processes, in the last few years. Figure 3 shows block diagram of a model predictive controller.
There are three main approaches to model predictive control, MAC (Model Algorithmic Control), which is based on system's impulse response, DMC (Dynamic Matrix Control), which uses the process step response samples, and GPC (Generalized Predictive Control), which is based on the process transfer function. In practice, it is easier to obtain step response samples rather than impulse response or a full transfer function, and therefore the DMC method is more popular. We use the DMC method in this research. The cost function is defined as:  is the reference input, the following filtered form is used as the tracking trajectory: The parameter α changes the place of the first order smoothing filter pole; the smaller α the faster output will become. It has been shown that system robustness can be decreased by the reduction of α and increment of the  [16]. Expanding the summations and substituting the quadratic forms, the cost function in (9) can be rearranged in a matrix form as There is no pure time delay in the model, therefore, 1 N is zero, then: In forming (14) which is called programmed MPC. If the future set points are unknown, d y is assumed to be constant during the prediction horizon, i.e.: which is called non-programmed MPC. For an LTI system, without any constraints on output or control signal, the above optimization problem has the following closed form:  (19) and i g s are the step response samples. The matrix + G is a Toeplitz matrix consisting the step response samples as shown in (19). There is no closed form solution for the formulated constrained optimization problem. Hence, an online optimization algorithm (a genetic algorithm) is applied to solve the problem, which will be discussed later. The model output is: where N is the number of system step response samples reaching to steady state or equivalent impulse response steps which lead to zero; and N g is the system DC gain.

THE MODIFIED DMC
If the system has any poles close to the origin, the step response will be very slow and the required N is very large. A system including integrator never reaches to the steady state (this case exists in the set of linearized models of the MMA reactor) and N approaches infinity. Hence, instability occurs. This is one of the limitations of the standard DMC formulation, making it only applicable to open loop stable system [17]. To overcome this, one practical solution is as follows. We have where Past Y is the effect of past inputs on the future system outputs without considering the effect of present and future inputs. Consequently, Past Y can be calculated by setting the future u ∆ s equal to zero and solving the model P steps ahead.
As seen in (15) and (19), + G and + ∆U are independent of N . The only thing determined by N is the dimension of − G , which is omitted now. Therefore, using this technique, DMC computations become independent of N . This modified formulation can be used for marginally stable and unstable plants alike.
As discussed in the previous section, we have used a piecewise linear model of the MMA polymerization process. In our application, since the valid model for future time steps may change through the course of predictions, in the computations of Past Y , in order to predict the future model outputs, the corresponding valid models are used. In other words, the model used for output prediction is scheduled through the prediction horizon, exploiting the developed trajectory linearized model. This virtual model switching is utilized to calculate the optimal control sequence in each time step. Apparently, another model switching is also applied through the course of time simulation of the closed-loop control system. Furthermore, since the whole desired temperature profile is known a priori, the programmed MPC is used, utilizing the known future desired temperature trajectory.

4 B GENETIC ALGORITHM OPTIMIZATION
In general, there is no closed form solution for MPC, except the linear time-invariant unconstrained case.
Otherwise, the optimization problem should be solved numerically. If the solution space is convex, sequential quadratic programming techniques (SQP) could be used. Otherwise, either the optimization problem should be a convexified through approximations and relaxation methods or a global optimization algorithm must be used.
Genetic algorithms (GA), are randomized global searching methods developed from the evolution rule in ecological world (the genetic mechanism of survival of the fittest). They have internal implicit parallelism and better optimization ability. By the optimization method of probability, they can automatically obtain and instruct the optimized searching space and adjust the searching direction autonomously [18]. Genetic algorithms, first introduced by Holland [19], are robust global random search methods. These methods are founded based on the Darwinian concept of natural selection and evolution [20], and have been used extensively in optimization and control [21][22][23].
Coding is essential for the GA optimization. Coding is a mapping from solution space to a finite length strings set. Binary coding is the most generally used method. In this method, each point of the solution space is coded as a binary string, which is a permutation of 0 and 1. Each 0 or 1 is called a gene and the string is called a chromosome. Every potential solution, , which is a point in the solution space, is coded as a binary word by length of , is a u N binary word, called a subchromosome. Figure 4 shows the binary coding of the + ∆U vector as a chromosome.
At the beginning, an initial population that consists of some potential solutions is randomly selected. The final solution is concluded by an iterative method that leads the initial population to an evolutionary one. In each iteration, the next population is generated by applying the crossover and mutation operators to the selected individuals. Their offspring makes the next population. Their selection probability depends on their fitness function that is a measure of goodness. The encoding method and fitness function definition are the only links between the physical problem and GA optimization. The following fitness function is applied: In this stage, pairs of the k-th population, k p , are selected to reproduce their offspring. The tournament selection strategy, which is a stochastic method, is applied to select each parent. In this method, two individuals are randomly selected and their best (according to their fitness value) is the winner of the tournament, then, the parents are returned to k p .

c. Crossover
This operator crosses parents to produce the new offspring by gene interchanging. Crossover may occurs in one or more positions in chromosomes. Researchers have suggested several crossover operators such as one point, multipoint and uniform crossover. In this study, a multipoint crossover operator is used to interchange genes between subchromosomes of parents. Figure 5 shows the multipoint crossover genetic operator. The crossover sites are determined randomly.

d. Mutation
In this phase, a random gene from chromosome is selected and its value will be changed. To do so, a random number between 0 and 1 is generated and compared with the predetermined mutation probability (  Fig.4. Binary coding of + ∆U as a chromosome directly transferred to the next one to prevent deterioration. The details of optimization algorithm are depicted in Figure 6. In the controller simulation, the optimization is conducted on a population of 50 chromosomes. As mentioned earlier, GA optimization belongs to the family of randomized search algorithms. It is worth mentioning that, in general, there are no theoretical proofs of the speed of convergence of randomized algorithms (convergence within certain time frame, if any). However, since the cost function is convex and the chemical processes are generally slow, allowing to have sampling periods in the order of tens of seconds, the speed of convergence should not be a concern here, provided that the population size and the evolution parameters are selected properly. At the same time, premature convergence is avoided by ensuring good genetic variation. The genetic variation can be regained by using a large enough population size and also by mutation [25]. On a separate note, in fast systems with millisecond sampling times, the optimization algorithm may not converge within the sampling interval and the optimization computation might be halted. In such cases, especially in mission-critical and safety-critical applications, hybrid algorithms must be utilized using FSQP (Feasible SQP) solvers, to ensure feasibility (satisfaction of all constraints) of the optimizer in each and every iteration even before convergence.

5 B SIMULATION RESULTS
The population is the foundation of evolution of a genetic algorithm. The character of the population decides the search capability of the genetic algorithm. The astringency of the genetic algorithm is determined by the astringency of the population [24].
First, the effects of population size and the number of generations on controller performance are studied. As it is shown in the Our several simulations show that the GA optimization algorithm converges well within the selected MPC sampling period in all runs. Figure 8 shows the simulation results of controller performance with a population size of 50 and the number of generations equal to 750. The top plot shows the systems output verses the desired thermal trajectory and exhibits controller's high performance. The middle plot is the control signal, which as seen in the plot, satisfies the constraint. The bottom plot is the tracking error, i.e. the difference between the actual and desired outputs. The DMC controller provides integral action, capable of rejecting step disturbance. The ability of controller to reject noise and disturbance is shown in Figure 9, in which a step output disturbance and a zero mean white Gaussian measurement noise with a variance of ) are applied. As seen, the controller has good disturbance rejection performance and the average absolute tracking error is less than . Compared to the previous results, the adaptive PI control in [13] has a 2° C average error and the Generalized Takagi-Sugeno-Kang fuzzy controller proposed in [26] has a 1°C average error; demonstrating the superior performance of the MPC.

CONCLUSION
A sequential piecewise linearized model based predictive controller based on the DMC algorithm was designed to control the temperature of a batch MMA polymerization reactor. Using the mechanistic model of the polymerization, a transfer function was derived to relate the reactor temperature to the power of the heaters. The coefficients of the transfer function were calculated along the selected temperature trajectory by sequential linearization. A genetic algorithm (GA) is applied for the cost function optimization DMC. The simulation result of controller performance shows that the tracking the profile, noise and disturbances rejection is very good.