CHEMICAL ENGINEERING TRANSACTIONS VOL. 70, 2018 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Guest Editors: Timothy G. Walmsley, Petar S. Varbanov, Rongxin Su, Jiří J. Klemeš Copyright © 2018, AIDIC Servizi S.r.l. ISBN 978-88-95608-67-9; ISSN 2283-9216 Multi-objective Optimisation of MTBE Reactive Distillation Process Parameters based on NSGA-II Wu Xiao, Yi Zhang, Xiaobin Jiang, Xiangcun Li, Xuemei Wu, Gaohong He* State Key Laboratory of Fine Chemicals, R&D Center of Membrane Science and Technology, Dalian University of Technology, Dalian 116024, Liaoning, China hgaohong@dlut.edu.cn Optimisation of reactive distillation (RD) process parameters is very important and complex because of a large number of variables and a strong coupling between them. A hybrid method of orthogonal numerical test and non-dominated sorting genetic algorithm II (NSGA-II) was presented for optimisation of the methyl tert-butyl ether (MTBE) RD process parameters such as the reflux ratio, the number of trays etc. on maximizing the conversion of isobutene and minimizing TAC of MTBE RD process simultaneously. Firstly, Aspen Plus was used to simulate MTBE RD process and the key parameters were selected by univariate sensitivity analysis. Secondly, the orthogonal test method was used to design the experiment scheme, and correlation formulations between the process parameters and objective functions (conversion of isobutene and total annual cost) were determined by second order regressions model according to the results of orthogonal numerical tests. Finally, based on the fitting function equations, NSGA-II was used to optimize the MTBE RD process parameters, and a non-dominated solution set was obtained and presented. The results indicate that the combination of orthogonal numerical test and the NSGA-II can not only reduce the number of experiments effectively, but also get the optimal and effective process parameters set for process design and operation. 1. Introduction Reactive distillation (RD) is a process that combines reaction and distillation into one step (Srinivas et al., 2010). Compared to the classical serial arrangement of unit operations, RD has the potential to decrease the dimensions of the equipment and to increase the degree of heat integration (Urselmann et al., 2011). During the past few ten decades, many innovative design methods for RD columns have been put forward (Chang and Lee, 2017). RD has been widely applied in the process of etherification, esterification and alkylation in the past (Estrada-Villagrana et al., 2006). The effect of equipment and operating parameters on RD process performance indicators (such as the conversion of key reactant and the composition at the bottom/top of column in RD) is more sensitive compared with traditional reaction or distillation process. For example, the variation of reflux ratio can result in not only the change of liquid composition of distillation tray, but also the variation of contact conditions between liquid and catalysts. Since the complex coupling interaction between reaction and distillation, the minor variations of process parameters such as feed position, theory plate number, reflux ratio, column pressure and feed temperature can bring strong variations of product purity and total annual cost (TAC). Therefore, the optimisation of RD process parameters is vital to industrial application, energy saving and consumption reduction (Alireza and Sirous, 2011). The simplifying treatment is needed during simulating RD process because it is complex. With the development of chemical process simulation software, the RD simulation based on strict thermodynamic equations can be achieved by Aspen Plus (Santoso et al., 2009). The RD process is highly nonlinear, and also contains integer variables, such as the plate number, feed location and other parameters, so the reactive distillation process optimisation is a mixed integer nonlinear programming (MINLP) problems. Noshadi et al. (2012) develop an optimal continuous process to produce fatty acid methyl esters (biodiesel) from waste cooking oil in a reactive distillation column catalysed by a heteropolyacid. Fatty acid methyl ester (FAME) yield was the response function, and Response Surface Methodology (RSM) was used to design the experiment and analysed four operating parameters: total feed flow, feed temperature, DOI: 10.3303/CET1870271 Please cite this article as: Xiao W., Zhang Y., Jiang X., Li X., Wu X., He G., 2018, Multi-objective optimisation of mtbe reactive distillation process parameters based on nsga-ii , Chemical Engineering Transactions, 70, 1621-1626 DOI:10.3303/CET1870271 1621 reboiler duty and methanol/oil ratio. Lu et al. (2017) Bat algorithm (BA) was used to optimize the RD for the production of methyl acetate (MeAc). Based on the link between Matlab and Aspen Plus. BA can find the global optimal solution within less computation time than other stochastic algorithms or sequential optimization. In recent years, the multi-objective optimization of reactive distillation has attracted widely attention. Non-Sorting Genetic Algorithm II (NSGA-II) is an effective algorithm for multi-objective optimization (Lv et al., 2017). NSGA- II and HYSYS software for thermodynamic calculation of RD column have been linked to optimize the effective parameters of RD column (Singh et al., 2005). The 9 parameters, such as feed ratio, reflux ratio etc., are selected as decision variables. The NSGA-II was employed for minimization of reboiler energy cost, maximization of n- butyl acetate molar flow as RD productivity, and maximization of methanol molar flow as non-reactive distillation column productivity. Multiple solutions set as optimal solutions are available for users. As a fuel additive, MTBE serve as a role of increasing octane number of gasoline. A large number of scientific works, which were developed for analysis and simulation of MTBE RD process, have been reported. Singh et al. (2005) developed the mathematical model of MTBE production by RD process. Santoso et al. (2009) used a Process Simulator to analysis MTBE Reactive Distillation Column. However, very few papers reported the optimisation of process parameters of MTBE RD. In this work, a study on MTBE RD process parameters optimisation was reported. The MTBE RD experiments were conducted with Aspen Plus. The parameters which affect the MTBE RD process were optimized through the combination of orthogonal numerical test and NSGA-II. The responses considered were the conversion of isobutene and the TAC. Pareto optimal solutions set, which was presented to satisfy the multi-choices for process design and operation, was obtained and reported. 2. Experimental work The numerical experiments were conducted to simulation of MTBE RD process by Aspen Plus. The selection of key decision parameters was performed by univariate sensitivity analysis. 2.1 Validation of MTBE RD model based on Aspen Plus To verify the modelling results of Aspen plus, the simulation output was compared with the corresponding results given in the literature for the same input parameters (Santoso et al., 2009). The results of MTBE RD model simulated with Aspen Plus under the above conditions show that the bottom temperature is 428.88 K and the purity of MTBE is 99.55%. The percentage errors of MTBE purity and bottom temperature are 0.31% and 0.54% respectively, which are within the acceptable range. Therefore, the simulation results are reliable. 2.2 Selection of key decision parameters Univariate sensitivity analysis is one of the key tools used to understand the effect of various parameters on the several performance indicators. In this study, the key decision parameters that influence the MTBE RD process are selected through univariate sensitivity analysis tools of Aspen Plus. According to the univariate sensitivity analysis, the key decision parameters for MTBE RD process are plate number (N), reflux ratio (R), feed position of methanol (SM), number of reactive trays (Nr), liquid holdup (L) and flow rate of methanol (FM). 2.3 Orthogonal numerical test The numerical experiments are conducted by using orthogonal design of experiments principle, which in turn reduces the number of experiments to be conducted in manufacturing organizations to obtain the results (Lee et al., 2003). Six factors and five levels orthogonal design L25(56) is used for the experimentation. The responses considered for the experiments are the total annual cost of MTBE RD process, and the conversion of isobutene. The value of each level for each factor is shown in Table 1. Table 1: The value of each level for each factor Factor level N R SM Nr L /m3 FM/kmol·h-1 1 2 3 4 5 15 16 17 18 19 4 5 6 7 8 5 6 7 8 9 5 6 7 8 9 3 3.2 3.4 3.6 3.8 620 640 660 680 700 The TAC for MTBE RD process (C) can be calculated by the Eq(1) (Seider et al., 2011). 1622 pUpHp CCCC  (1) Where pC = equipment cost of MTBE RD column, pHC = equipment cost of condenser and reboiler, pUC = utility cost. The parameters for calculation of capital cost and operating cost are given in Table 2. Table 2: The parameters for calculation of capital cost and operating cost item unit price/$·ton-1 Notes carbon steel 3230.00 DIN17155 cooling water 0.015 Inlet temperature: 298.15K Exit Temperature: 313.15K heating steam 8.20 437.45 K Table 3: Results of orthogonal numerical test Item N R SM Nr L FM TAC (106$/y) Conversion of isobutene 1 1 1 1 1 1 1 12.664 0.5295 2 1 2 2 2 2 2 12.410 0.7539 3 1 3 3 3 3 3 13.149 0.8485 4 1 4 4 4 4 4 13.755 0.9135 5 1 5 5 5 5 5 14.698 0.9587 6 2 1 1 3 4 5 10.761 0.8205 7 2 2 2 4 5 1 11.508 0.8315 8 2 3 3 5 1 2 12.295 0.8908 9 2 4 4 1 2 3 15.552 0.8017 10 2 5 5 2 3 4 16.001 0.8746 11 3 1 4 5 2 4 9.324 0.9315 12 3 2 5 1 3 5 9.612 0.8691 13 3 3 1 2 4 1 14.827 0.7001 14 3 4 2 3 5 2 14.900 0.8220 15 3 5 3 4 1 3 15.390 0.8937 16 4 1 5 2 5 3 10.234 0.8608 17 4 2 1 3 1 4 12.310 0.8170 18 4 3 2 4 2 5 12.851 0.9116 19 4 4 3 5 3 1 14.046 0.8681 20 4 5 4 1 4 2 17.146 0.7767 21 5 1 4 4 3 2 9.938 0.8692 22 5 2 5 5 4 3 11.164 0.9192 23 5 3 1 1 5 4 16.730 0.6235 24 5 4 2 2 1 5 15.804 0.8183 25 5 5 3 3 2 1 16.358 0.8247 The conversion of isobutene can be calculated by using the Eq(2) )(1 xf = 1- (FTOP,iso + FBOTTOM,iso)/FFEED,iso (2) Where )(1 xf is the conversion of isobutene; FTOP,iso is the flow rate of isobutene at the top of column, kmol·h-1; FFEED,iso is the feed flow rate of isobutene, kmol·h-1; FBOTTOM,iso is the flow rate of isobutene at the bottom of column, kmol·h-1. Results of orthogonal numerical test at the constraint that the purity of MTBE at the column bottom is greater than 97 % are given in Table 3. 3. Mathematical model The conversion of isobutene has an important influence of the TAC of MTBE RD process. The conversion of isobutene and TAC are selected as the key objectives. The polynomial second order is developed to establish the relation between the process parameters and response of the conversion of isobutene in MTBE RD process is described as follows: 1623    ji i j i k i k i i xxxxy ij 2 1 ii 1 i0 ββββ (3) Where β is coefficient, used in the above model can be calculated by Eq(4): YXXX TT 1)(β  (4) Where TX is the transpose of matrix X, and Y is the matrix of measured response under study. For calculating the coefficients, least square method is used, and a computer program is compiled. The final mathematical model formation developed the conversion of isobutene is given as: 2 665 2 56454 2 4 635343 2 3625242 32 2 26151413121 2 16543211 0007.0001198.09105.0000094.024704.006837.0 00217.021222.03926.000715.000102.0009577.05238.00 4526.002978.000025.000204.306378.000431.000194.00 3932.004573.008094.121233.004292.81.0143600714.6090857.0)( xxxxxxxxx xxxxxxxxxxxxx xxxxxxxxxxxxx xxxxxxxxf     (5) The linear model is developed for the TAC is given based on the same method: 653212 00519.042921.047275.050249.114212.091501.9)( xxxxxxf  (6) Where 1x number of trays N, 2x reflux ratio R, 3x feed position of methanol SM, 4x number of reactive trays Nr, 5x liquid hold up L, 6x flow rate of methanol FM. The model adequacy was checked by using the coefficient of correlation value. The coefficients of correlation are calculated for both models and are found to be 0.9994 for the conversion of isobutene and 0.9327 for the TAC, which shows the high correlation that exist between the models and experimental results. 4. Optimisation of MTBE RD process parameters with NSGA-II 4.1 multi-objective optimisation problem In optimisation of MTBE RD process parameters, NSGA-II is used. Two-objective optimisation is carried out in this study. The objectives set for the present study is as follows: (1) Maximization of the conversion of isobutene ( )(1 xf ). (2) Minimization of TAC of MTBE RD process ( )(2 xf ). The optimisation problem is constrained by the process-specific constraints, which are given as the upper and the lower bounds on the conversion of isobutene. The lower bound is zero. The upper bound of 0.95 is defined in order to consider the fact that the equation )(1 xf cannot strictly represent the relation between the conversion of isobutene and key decision parameters when the conversion of isobutene is greater than 0.95. The TAC is greater than zero. This constraints of the multi-objective problem can then be expressed as follows: 0< )(1 xf < 0.95; 0< )(2 xf ;15 ≤ 1x ≤ 19; ;1 Nx  4 ≤ 2x ≤ 8; 5 ≤ 3x ≤ 9; ;3 Nx  5 ≤ 4x ≤ 9; ;4 Nx  3 ≤ 5x ≤ 3.8; 620≤ 6x ≤700. 4.2 NSGA II algorithm NSGA II algorithm is used to solve the multi-objective problem. There are two key concepts in NSGA-II: a fast non-dominated sorting of the population and a crowding distance. The solution procedure of NSGA-II is described as follows: (1) Initialize population (pop), number of generations (gen), number of decision variables (V), and number of objectives (M). (2) Evaluate objective functions. (3) Assign rank to each individual in population on the basis of non-dominance. (4) Sort the population using non-domination-sort. (5) Crowding distance calculation for each individual in the population. (6) For each generation: (a) Select the parents. Parents are selected for reproduction to generate offspring. The NSGA-II uses a binary tournament selection based on the crowded-comparision operator. 1624 (b) Generate offspring population by performing crossover and mutation based on the crossover and mutation probability. (c) Intermediate population. Intermediate population is the combined population of parents and offspring of the current generation. The population size is almost one and half times the initial population. (d) Non-domination-sort of intermediate population. The intermediate population is sorted again based on non- domination sort before the replacement operator is performed on the intermediate population. (e) Perform Selection. Once the intermediate population is sorted only the best solution is selected based on it rank and crowding distance. Each front is filled in ascending order until the addition of population size is reached. The last front is included in the population based on the individuals with least crowding distance. (f) Replace individuals in the population. Define that the number of tournament candidates is 2, the crossover distribution index and mutation distribution index both equal to 20, the crossover probability is 0.9, the mutation probability is 0.1. 5. Results and discussion The NSGA-II is programmed with Matlab and run successfully in a personal computer. The results of this multi- objective optimisation problem indicate that the objectives for the considered case are conflicting and this is hence presented as a Pareto-optimal front. According the calculation results of Figure 1 and 2, the generation of the NSGA-II is set as 500, and the population is set as 100. -0.95 -0.94 -0.93 -0.92 -0.91 -0.90 -0.89 -0.88 8.60 8.65 8.70 8.75 8.80 8.85 Population 100 and generation 100 Population 100 ang generation 500 Population 100 and generation 1000 A n n u a l to ta l co st (1 0 ^ 6 $ /y ) Inverse number of the conversion for isobutene -0.96 -0.94 -0.92 -0.90 -0.88 8.55 8.60 8.65 8.70 8.75 8.80 8.85 Population 20 and generation 500 Population 50 and generation 500 Population 100 and generation 500 A n n u a l to ta l c o st ( 1 0 ^ 6 $ /y ) Inverse number of the conversion for isobutene Figure 1: Optimization results for various generations Figure 2: Optimization results for various populations It can be seen from Figure 3 that the TAC decreases sharply with the conversion of isobutene in CD segment but smoothly in AC segment, and the Pareto solutions distributed within BC segment are preferable when both optimized objectives have to be satisfied simultaneously. The process engineer should select the proper solution according to his requirement. If the conversion of isobutene is more important than TAC, the process parameters will be selected from AB segment. On the contrary, if the engineer prefers the lower TAC, the parameters will be selected from CD segment. The process parameters should be selected from the BC segment when two objectives are equally important. -0.95 -0.94 -0.93 -0.92 -0.91 -0.90 -0.89 -0.88 8.60 8.62 8.64 8.66 8.68 8.70 8.72 8.74 8.76 Population 100 and generation 500 A n n u a l to ta l co st (1 0 ^ 6 $ /y ) Inverse number of the conversion for isobutene A B B C D -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 8 10 12 14 16 18 process parameters at the table 4 optimal process parameters with population 100 and generation 500 A n n u a l to ta l c o st ( 1 0 ^ 6 $ /y ) Inverse number of the conversion for isobutene Figure 3: Optimisation results for population 100 and generation 500 Figure 4. Comparison results of optimal solutions with the parameters in Table 3 1625 The optimisation output is compared with the corresponding results given in Table 3. The comparison result is shown in Figure 4. It can be seen from Figure 4 that the optimal MTBE RD process parameters with the results of population 100 and generation 500 owns higher conversion of isobutene at the precondition of lower TAC compared with the process parameters at the Table 3. This indicates that the combination of NSGA-II and orthogonal numerical test not is effective for getting Pareto frontier in multi-objective optimisation of MTBE RD process parameters. 6. Conclusions A hybrid approach combining orthogonal numerical test and the NSGA-II is proposed for optimisation of RD process parameters in the work. The orthogonal numerical test is used to design the experiment schemes. Then the polynomial second order is developed to establish two relation equations between the parameters and two different objectives. Finally, the optimal MTBE RD process parameters are obtained through NSGA-II based on multi-objective optimisation strategy. A case study shows that, based on the new method, higher conversion of isobutene and lower TAC can be achieved in the optimal process parameters. More importantly, the proposed design method of combination of orthogonal numerical test and the NSGA-II not only can reduce the number of experiments effectively, but also can get the optimal and effective process parameters set for process design and operation. The effectiveness of application can be further improved by using other methodologies to select key decision parameters such as Plackett-Burman design. Acknowledgments We grateful thank the financial support from National Natural Science Foundation of China (21676043); Program for Changjiang Scholars (T2012049); Fundamental Research Funds for the Central Universities (DUT17JC33, DUT17ZD203, DUT16TD19) ; MOST innovation team in key area (No. 2016RA4053), Education Department of the Liaoning Province of China (LT2015007) . References Alireza B., Sirous S., 2011, Multiobjective optimisation of reactive distillation with thermal coupling using non- dominated sorting genetic algorithm-II, Journal of Natural Gas Science and Engineering, 3(2), 365–374. Chang H.N.,Lee H.Y., 2017, Innovative Design of Diphenyl Carbonate Process via One Reactive Distillation with a Feed-Splitting Arrangement, Chemical Engineering Transactions, 61, 445-450. Estrada-Villagrana A.D., Quiroz-Sosa G.B., Jiménez-Alarcón M.L., Alemán-Vázquez L.O., Cano-Domínguez J.L., 2006, Comparison between a conventional process and reactive distillation for naphtha hydrodesulfurization. Chemical Engineering and Processing Intensification, 45(12), 1036–1040. Lee, K.Z., Chuang W.C., Ho S.Y., 2003, A non-parametric image segmentation algorithm using an orthogonal experimental design based hill-climbing. Intelligent Data Engineering and Automated Learning, 2690, 1076– 1081. Lu J.,Tang J.,Chen X.,Cui M.,Fei Z.,Zhang Z.,Qiao X., 2017, Global Optimization of Reactive Distillation Processes Using Bat Algorithm. Chemical Engineering Transactions, 61, 1279-1284. Lv J.F., Jiang X.B., He G.H., Xiao W., Li S., Sengupta D., Mahmoud M.E., 2017, Economic and system reliability optimization of heat exchanger networks using NSGA-II algorithm, Applied Thermal Engineering, 124, 716– 724. Noshadi I., Amin N.A.S., Parnas R.S., 2012, Continuous production of biodiesel from waste cooking oil in a reactive distillation column catalyzed by solid heteropolyacid: Optimisation using response surface methodology (RSM), Fuel, 94, 156–164. Santoso H., Bao J., Lee P.L., 2009, Operability Analysis of MTBE Reactive Distillation Column using a Process Simulator. Chemical Product and Process Modeling, 4(3), 1934-2659. Seider W.D., Seader J.D., Lewin D.R. (Ed), 2004, Product and process design principles: synthesis, analysis, and evaluation. Wiley, New York, USA. Srinivas, S., Mahajani S.M., Malik. R.K., 2010, Reactive Distillation for Fischer-Tropsch Synthesis: Simulation- Based Design Methodology Using Aspen Plus, Industrial & Engineering Chemistry Research, 49(20), 9673– 9692. Singh, B.P., Singh R., Kumar M.V.P., Kaistha N., 2005, Steady-state analyses for reactive distillation control: An MTBE case study, Journal of Loss Prevention in the Process Industries, 18(4-6), 283–292. Urselmann M., Barkmann S., Sand G., Engell S., 2011, Optimisation-based design of reactive distillation columns using a memetic algorithm, Computers and Chemical Engineering, 35(5), 787–805. 1626