Reactor Network Synthesis Using Coupled Genetic Algorithm with the Quasi-linear Programming Method

This research is an attempt to develop a new procedure for the synthesis of reactor networks (RNs) using a genetic algorithm (GA) coupled with the quasi-linear programming (LP) method. The GA is used to produce structural configuration, whereas continuous variables are handled using a quasi-LP formulation for finding the best objective function. Quasi-LP consists of LP together with a search loop to find the best reactor conversions (xi), as well as split and recycle ratios (yi). Quasi-LP replaces the nonlinear programming (NLP) problem, and is easier to solve. To prevent complexity and ensure an optimum solution, two types of ideal reactors, namely plug flow reactor (PFR) and continuous stirred tank reactor (CSTR), were considered in the network. Since PFRs require the introduction of differential equations into the problem formulation, a CSTR cascade was used instead in order to eliminate differential equations. To demonstrate the effectiveness of the proposed method, three reactor-network synthesis case studies are presented.


Introduction
The main objective of chemical processes is to convert raw materials into valuable products.Since, according to Smith 1 , the reactor network (RN) is the living heart of any industrial process, the proper design of this part results not only in the reduction of production expenses, but will reduce overall costs by affecting separation and utility costs.RN synthesis can be formulated as an optimization problem.In this problem, the decision variables are reactor volume, reactor type and configuration, yield, and selectivity resulting from various parallel and series reactions.Because of high nonlinearity of equations describing the system and complexity of kinetic equations, this problem is a hard optimization problem.This may explain why much fewer papers on this area of design are published than on other areas such as heat exchanger or separation networks 2,3 .
There are three different methods of design of RN in the literature.The first is used in industrial processes and applies to modifications of existing plants, it is highly dependent on the intuition and experience of the designer, and is based on some heuristics developed mainly on the basis of the work of Douglas 4 .
The second method is based on mathematical programming, which seeks to optimize configuration.This method is divided into two categories.The first one optimizes the structure 2,5,6 , while the second gives a systematic way to reach an optimum configuration 7,8 .In the structure optimization, a primary super-structure is initially proposed and an optimal sub-network that optimizes a desired variable is derived from the initial network.However, the primary structure may not contain all potential sub-structures 2,3,9 .
The third method of design is based on the systematic generation of the process flow sheet or part of it.Here, the objective is to find all possible concentration domains resulting from various reactor systems, and to consider all defined reactions 10,11 .Therefore, a functional representation is used to model all reactions and mixing.This is based on concentration state space and is called the attainable region (AR) 12,13 .This method was first devised by Horn 10 .After Glasser et al. 11 had defined a systematic graphical method, this method became popular in the field of design of reactor networks.The main characteristic of the two latter methods is that they can be considered as mathematical programming methods, and the problem can be solved to find the optimum process flow diagrams.
The attainable region method, as an applicable method, was used in the works of Glasser et al. 11 and Glasser and Hildebrandt 14 .In their method, the reactors were assumed isothermal, while change in volume inside the reactors and in mixing was ignored.
McGregor 12 outlined a trial and error search procedure in state space to find a candidate attainable region.He considered boundaries for the region, and the constraints were examined for the region.If any constraint deviated, another candidate region was examined until all constraints were satisfied.Although the method did not require considering reactor configuration or intuition of the designer, it was applicable only to small, two-dimensional problems.The difficulties stemmed from the fact that imagination of geometries requiring more than three parameters is difficult if not impossible, and that this method was time-consuming.Rooney et al. 15 proposed a method to overcome the first difficulty, which was a four-stage algorithm for finding the attainable region.The region was found by producing a set of two-dimensional attainable regions under steady sub-spaces.These sub-spaces could be produced by special interpretation of CSTR and PFR reactors.Mathematical programming was also used in this method.To solve the second difficulty, i.e., time consumption, the research was centered on computer structures.The result was Infinite DimEn-sionAl State-space (IDEAS), which appeared in the work of Burri et al. 16 In this method, the point in the boundaries of the region and in the concentration space could be found by Infinite -dimensional Linear Program (ILP).Manousiouthakis et al. 17 extended the Shrink-Wrap Algorithm (SWA) to state necessary and sufficient conditions for a point to be in the attainable region.Another extension of this method was published in the works of Zhou and Manousiouthakis, such that their methods could solve non-isothermal problems 18,19,20 .
At the same time, others tried to develop mathematical programming methods.Achenie and Biegler 21 presented an approach to convert the synthesis problem formulation into an optimal control formulation, solved using a gradient-based algorithm that employs successive quadratic programming and adjoint variables.Kokossis and Floudas 2 proposed a systematic method to solve the problem based on hyperstructures of reactors.They used a CSTR cascade instead of PFR.In this way, they prevented the appearance of differential equations in the formulation because of PFR reactors.They formulated the problem as a Mixed Integer Nonlinear Programming (MINLP).Balakrishna and Biegler 7 developed a method based on optimizing streams between various reaction media, which acted as a dynamic optimization problem.Bikic and Glavic 22 suggested an algorithm for producing reactor networks, which considered mixing and reaction.Cordero et al. 23 combined NLP with simulated annealing to design and optimize reactor networks by suggesting potentially feasible networks.Pahor et al. 24 presented a superstructure-based MINLP approach to the RNs in an equation-oriented environment, comprising isothermal and non-isothermal reaction problems.Revollar et al. 25 suggested a method for process integration, which considered control and economic approaches simultaneously.The problem was formulated as MINLP and a genetic algorithm was used to solve it.
Silva et al. 3 presented a calculation procedure to design isothermal RN.A superstructure was proposed with ideal CSTRs and PFRs with various configurations.The problem was formulated as a constrained NLP.The model was solved using IMSL/Fortran.Genetic Algorithm was used to define and configure reactors by maximizing yield or selectivity and minimizing reactor volume.
Labrador et al. 26 presented a novel framework for the optimization and synthesis of complex reactor networks that combined superstructure-based optimization, semantic models, and analytical tools.Pontes and Pinto 27 presented a MINLP model to synthesize Fenton reactor networks for phenol degradation.In this method, the network superstructures with multiple Fenton reactors were optimized with the objective of minimizing the sum of capital, operating, and depreciation costs of the effluent treatment system.For finding an optimal flow sheet structure and a robustly stable operating point that minimizes an overall cost function, Zhao and Marquardt 28 presented a nonlinear dynamic model for reactor networks synthesis with guaranteed robust stability.Their work extends the normal vector approach, which so far assumes a fixed flow sheet, and thus considers only continuous degrees of freedom for optimization.Hence, they generalize the method to integer degrees of freedom to support decision-making on the flow sheet level.Jin et al. 9 proposed a two-level procedure.In the inner level, reactor volume and flow rates were calculated using an LP.Outlet concentrations which were needed for this level were obtained by optimized search in the concentration space in the outer level.This search was done by combining the simulated annealing (SA) method and the particle swarm optimization (PSO) method.
This paper addresses the problem of the synthesis of RNs using a combination of genetic algorithm and quasi-linear programming (LP) method, in which GA, with the help of its operators (i.e.reproduction, crossover, and mutation) produces dif-ferent networks during the optimization process and finds the best structural modifications.The quasi-LP, which was formulated to find the best objective function value, consists of an LP and a search loop to find the best reactor conversions split and recycles ratios.The use of stochastic methods is justified by the power of these methods in handling the integer variables that appear in dealing with configurations  . Unti now, various forms of genetic algorithm (GA) 3 , simulated annealing (SA) 23,30 , Tabu Search (TS) 29,31 , and a combination of SA and PSO 9 together with optimization methods and attainable region, have been used to design RNs.The method used in this paper differs from the others in that it formulates the problem as a quasi-LP, and the configuration of the reactors is controlled by an address matrix which is compatible with GA operators.

Problem statement
Reactor networks synthesis is an optimization problem, therefore its solution requires input data.The problem is defined as follows: Given:

Methodology
The overall algorithm of the proposed method is shown in Fig. 1, in which the main objective is the maximization or minimization of the objective function (for example, minimum total reactor volume or total cost and/or maximum selectivity or the concentration of the desired product).In this figure, optimization is started with an initial population produced by the GA.In the next step, the networks are evaluated by the quasi-LP.This section consists of an LP and a search loop.To find optimal values of all variables by quasi-LP, these variables are classified into two groups.The first group consists of variables, such as concentrations and flow rates.The second group contains variables that make the equations nonlinear, such as conversion of each CSTR (x i ), the split ratios as well as recycle ratios (y i ).With known values of the second group variables in search loop, the first group variables can be obtained either by solving linear set of equations or by solving an LP, which is explained later.To start optimizing an objective function for a configuration produced by GA, some initial values are given to the second group of variables, then as explained above, the first group of variables is evaluated.In this way, a value can be found for objective function having values of all variables.Then the steps are repeated by modifying the values of the second group in the direction of optimizing objective function.
In this work, the conversion of each CSTR reactor, and split and recycle ratios are bound to [0.001 -0.999].It has been observed that whatever the lower amount of conversion is considered for each CSTR reactor in series, the results obtained from replacement of a PFR with a CSTRs cascade are better 2 .Therefore, the conversion of CSTR reactors in series is bound to [0.001 -0.3].
After determining the fitness of the first generation, the second generation is produced using genetic algorithm operators, which are reproduction, crossover, and mutation.The overall algorithm may be repeated for a specified number of generations or stopped when the desired solution is obtained.In this article, the former case is considered.
The following section explains the method for RN synthesis.

Representation of a network
In the present method, an RN is treated as a chromosome that has three genes, and in each gene, a maximum of two reactors is addressed.The number of parallel reactors in each gene as well as the number of genes in each chromosome can be increased easily, but this affects the number of governing equations and thus the runtime of the program.To show the location of reactors and their types in the network, a reactor address matrix (RAM) is defined, in which each gene comprises two rows of the matrix.Fig. 2 shows a network and its RAM with three genes and two reactors.In the RAM, the first column indicates split (2 for split and 1 nonsplit).The second column indicates reactor (2 for CSTR, 3 for PFR, 0 for no flow, and 1 for a stream).The third column shows recycle stream for PFR reactor.Note that the recycle stream is defined only for PFR.If necessary, a bypass stream can be considered.In the figures, C and P stand for CSTR and PFR, respectively.Furthermore, to eliminate differential equations in the model, the PFR is replaced by a CSTR cascade, which is shown in Fig. 3.It should be noted that a bypass stream is considered in any structure for managing environmental and safety precautions, therefore a bypass stream is not considered in the reactor address matrix (RAM).This representation is especially suited for GA operators, and always creates feasible structures.Three operators are considered to obtain the best objective function: reproduction, crossover, and mutation.

GA operators
The GA operators are reproduction, crossover, and mutation.
(a) Reproduction: In each population, some structures are reproduced by roulette wheel according to their fitness.These structures are exactly copied to the next generation according to the survival rate, which is set to 40-50 %.For the robustness of the algorithm, elitism is added to copy the best solution to the next generation.The number of elites depends on the size of the population but ranges between 1 and 5.
The number of networks in the initial population depends on the size of the problem.Note that the initial population is produced randomly.It is clear that the number of generations is proportional to the size of the initial population and the chromosome length.So, in this method, normally 40-60 generations are employed to reach the stopping criteria of the algorithm.
(b) Crossover: To produce a new generation by crossover operator, both single-point and two-point crossovers are employed.After selecting two networks by roulette wheel as parents, a random number is generated to decompose the chromosome into several pseudo-parents (four parents for single-point crossover and six for two-point).These parents are then combined together to produce an offspring.In this paper, the rate of crossover is considered to be 50-60 %.
(c) Mutation: The definition of the mutation and its rate plays an important role in the convergence of the algorithm.Actually, an inappropriate definition of this operator leads to premature convergence of the GA.Therefore, mutation should be defined in a manner that the GA has sufficient diversity during its evolutionary procedure.In this method, a random number is generated for each chromosome in the network, and this number determines if that chromosome should be mutated.If mutated, a part of the structure of the chromosome is changed.A variable rate of mutation for each generation, which had been used by Ravagnani et al. 32  Mass and mole balance equations around mixing and split points, as well as reactors will be as follows, according to Figs. 2 and 3. Mass balance for point (1) as a split point (constant density): for point (1-1) as a split point: for point (1-2) as a mixing point: for the reactor in the first gene: for points (2) and ( 3): for point (3-1) as a mixing point: for point (3-2) as a split point: for PFR or CSTRs in series in the third gene: for point (4) as a mixing point: Mole balance for component A for point (1): for point (1-1): for point (1-2): for CSTR in the first gene: for points 2 and 3: for point (3-1): for point (3-2): for PFR or CSTRs in series in the third gene: (18) ...
(1 ) for point (4): Mole balance for component B (representing other species other than A) for point (1): for point (1-1): for point (1-2): for CSTR in the first gene: for points (2) and (3): for point (3-1): for point (2-3): for PFR or CSTRs in series in the third gene: (27) ) for point (4): for the CSTR volume: for the PFR volume: Here, PFR has been replaced by a CSTR cascade.Therefore, volumes of CSTR are calculated to obtain PFR volume, as follows: (Fig. 3) In the above equations ν ij , c Aij , c Bij , x Ai , y i express j th stream flow rate, concentration of A, concentration of B, conversion of A, and split ratio of i th gene, respectively, and x repre- sent j th stream flow rate, concentration of A, concentration of B, and conversion of A, respectively, in the i th gene at k th CSTR reactor in series in place of PFR (Fig. 3).Also V C and V P are CSTR and PFR volumes, respectively, and V C k is the k th CSTR volume representing PFR.Rate of reaction is r A , and F represents reaction function which will differ according to the reaction equation.

Quasi-LP procedure
Considering the governing equations of the reactor network in Fig. 2, (Eqs. 2 to 31), the quasi-LP formulation includes the following items: Firstly, in the outer loop, initial values are assumed for continuous variables of the second group.With these variables, continuous variables of the first group are calculated according to the following steps: Step one -Finding the volume flow rates of streams and concentration of A inside them: Ignoring volume change by mixing or passing through reactor and having estimated the values of split ratios and recycle ratios, as well as values of conversions of CSTR obtained from the search loop, the equation of flow rates and concentration of A in streams (i.e.Eqs. 2 to 19) will be linear.Therefore, by calling an LP procedure the objective of which is maximizing conversion of species A in RN, is defined as: All values of flow rates and concentrations of species A in the streams of RN are obtained.The reason for using this objective function is to minimize the amount of unreacted materials in the output of the reactor network.This decreases the cost of separation of components, as well as the recycle stream of unreacted materials to the reactor (or reactors).Furthermore, this step has been programmed to embed any constraint (for any reason, such as safety or environmental considerations) on conversion of A and flow rate in each stream.
Step two -Obtaining other materials concentrations: Equations 20 to 28 show that other concentrations depend on A and volume flow rates, which have been calculated in step one.In addition, de-pending on the type of reactions, some concentrations may depend on other concentrations.For example, for the Van der Vusse reaction: (33) the outlet concentration of B, C and D in any CSTR could be expressed as: where in and out correspond to inlet and outlet of CSTR, respectively.Therefore, to calculate the concentrations of other species in RN, a set of nonlinear equations will be treated (see Equations 34 to 36 for Van der Vusse reaction).Solving simultaneously the set of nonlinear equations for all concentrations of other species in RN may be time-consuming if not impossible.
However, regarding the encoding of RN by genes, which is consequential, it is possible to solve the set of nonlinear equations by Newton's method for variables of each gene.Having obtained the values of the gen, the values are used as an initial guess for the input variables of the next gene, and so on.The ability of Newton's method to obtain values of the concentrations of other species is reliable, as demonstrated in the case studies.
In this step, it is possible to consider any constraint on concentration of materials.After solving the nonlinear equations and obtaining the values of variables if there are constraints that are not satisfied, a penalty is considered, which is added to the value of the objective function evaluated in step three.
Step three -Calculation of the overall objective function: The objective function is formulated based on the conditions of the problem.This function may be based on maximum concentration of an intermediate or maximum selectivity or minimum volume of reactors for concentration of the desired product, etc.If there are any penalties in the previous steps, they are added to the function here in this step.
After Step 3, new values are considered for conversion in each reactor split and recycle ratios, and all three steps are repeated until a new objective function is obtained, and depending on the value of the new objective function, the new values of the conversions as well as split and recycle ratios are determined.These operations continue until satisfactory values are obtained for the optimum objective function.

Illustrative examples
In order to validate the proposed approach,three complex reaction schemes from the literature were considered.The first example is the Trambouze reaction, the second is the Van der Vusse reaction, and the third is conversion of α-pinene reaction.
The original codes were compiled to C language, and the final codes were loaded on a computer center with five parallel computers of 4 GB of ram memory and 3.4 GHz processor with a master processor.The optimization time that corresponds to each case study is about 50 minutes, which may be reduced by increasing the number of slave computers and improving their characteristics.
The number of CSTRs in series to simulate PFR was 50.Although increasing this number increases precision, the time and volume of calculations increase dramatically, therefore, this number seemed reasonable.The following assumptions were also made: 1. Reactions take place in liquid phase 2. Stream volumes remain constant in mixing 3. Reactors are isothermal.

Example 1 -Trambouze Reaction Scheme
The Trambouze reaction scheme is defined by the following reactions, which involve four species, and is composed of three parallel reactions where C is the desired product. (37) The data of the reaction is shown in Table 1.The objective is to maximize the selectivity of component C to component A, which is defined as: where x C and x A are mole fractions of C and A, respectively.Fig. 4 shows the best reactor network obtained by the proposed method, which includes a CSTR and a bypass that is in series with a plug reactor with recycle.The maximum selectivity is 0.5.One of the advantages of GA is that it does not produce only one solution.The second and third best RN is shown Fig. 5.
Table 2 compares optimization results with those in the literature.As shown in this table, the value for the objective function is the same; however, the total volume of the reactors obtained by this method compared to those in the literature is much lower.The Van der Vusse reaction scheme is defined by the following reactions, where B is the desired product.

Ta b l e 1 -Information of Trambouze reaction -case study one
(39) The reaction data is shown in Table 3.The objective is to maximize the outlet concentration of B.
Fig. 6 shows the optimum RN obtained by the proposed method, which consists of two parallel CSTR in series with a PFR.The maximum outlet concentration of B is 0.68754 mol L -1 .8, which contains a CSTR with bypass in series with a PFR.For this structure, the maximum outlet of B is 0.68752 mol L -1 .
Table 4 compares these results with those in the literature.As shown in this table, both solutions gave better results as regards B concentration and minimum reactor volume.

Example 3 -Conversion of α-pinene Reaction Scheme
This scheme is a complex reaction first introduced by Fuguitt and Hawkins 33,34 .Stewart and Sorensen 35 performed a research on the reaction constants and obtained the constants at 580 °C.The reaction is as follows 21 : (40) Table 6 compares the results with those in the literature.As can be observed, the selectivity compared to the Ta b l e 3 -Information of Van der Vusse Reaction-case study two

Conclusion
In this paper, a new effective method for synthesizing reactor networks is proposed based on GA-quasi-LP method.To show the location of reac-tors and their types in the network, a reactor address matrix is defined which can be coded with GA operators.In this RAM, each gene reserves two rows of the matrix, so RAM can represent nearly all structures.To reduce complexity, only two types of reactors, namely CSTR and PFR, were used.To remove differential terms in material balance equations for PFR, the PFR is replaced with a CSTR cascade.Although the error is reduced by using more CSTRs in series, the volume of calculations also increased.In this study, 50 CSTRs in series were used for simulating a PFR, which seemed reasonable.
Instead of obtaining all continuous variables simultaneously inside the quasi-LP, some variables that cause the equations to become nonlinear (such as split and recycle ratios, as well as conversions) were estimated in a search loop, while other variables (such as flow rate and concentration of component in stream) were obtained in the LP using mass and mole balance equations.Using a threestep procedure inside the LP section, it is possible to consider any constraint in addition to obtaining optimum values of concentrations and flow rates.In Ta b l e 5 -Information of α-Pinene Reaction-case study three v(PureA)  this way, the synthesis of the reactor network is con verted from complex MINLP model to GA-Quasi-LP model, which is much easier to solve, as well as increases the probability of reaching a global optimum.
The new method can be applied easily to non-isothermal and industrial cases, and this is the subject of future research work.The benefit of this kind of representation is that only feasible networks are produced during the optimization process.
One of the advantages of this method is that it needs no initialization because of the use of the GA in structural optimization and simplex method in the LP formulation.The only variables that require an initial guess are the conversion of component Ain each CSTR (x i ), and split and recycle ratios (y i ).In this paper, an initial value of 0.5 with upper and lower boundary of 0.001 and 0.999 was considered for all split and recycle ratios, as well as CSTRs conversions.In addition, a value of 0.1 with upper and lower boundary of 0.001 and 0.3 was used for conversion of CSTRs in series representing a PFR, since better results are obtained for lower conversions 2 .

F i g . 1 -
Overall algorithm for synthesis of RN: (a) GA and (b) Quasi-LP procedure for synthesis of RN , was used.When minimizing the objective function: Mu = Mu min + (Mu max -Mu min ) exp [-10(OB w -OB b )/OB w ] (1a) When maximizing the objective function: Mu = Mu min + (Mu max -Mu min ) exp [-10(OB b -OB w )/OB b ] (1b) In Eq. (1a and b), Mu is the mutation rate, calculated in each generation.Mu min and Mu max are the minimum and maximum mutation rates allowed, and their values are 10 % and 100 %, respectively.OB w is the value of the worst objective function, and OB b is the value of the best one in the previous population.Mathematical model of superstructure and LP procedure Mathematical model of superstructure

F
i g . 2 -RN with three genes and two reactors, and its RAM F i g .3 -Simulation of PFR by CSTRs in series-(a) PFR and (b) CSTRs in series

F i g . 4 -
Optimum network obtained for the first case study F i g .5 -Other optimum networks obtained for the first case study

F i g . 6 -Fig. 7
Fig.7shows the maximum and average mole concentrations of populations (for component B) vs. iterations.As shown in this figure, this structure was achieved after 60 iterations.The second best solution is shown in Fig.8, which contains a CSTR with bypass in series with a PFR.For this structure, the maximum outlet of B is 0.68752 mol L -1 .Table4compares these results with those in the literature.As shown in this table, both solutions gave better results as regards B concentration and minimum reactor volume.
the α-pinene (C 10 H 16 ), B is the α-and β-pyronene (C 10 H 16 ), C is the dimmer (C 20 H 32 ), D is the allo-ocimene (C 10 H 16 ), and E is the dipentene (C 10 H 16 ).The reaction data is shown in Table5.The objective function is the selectivity of C over D, defined as: molar fraction of C/molar fraction of D, such that outlet concentration of D from RN is at least 0.01 mol L -1 .Fig.9shows the optimum network obtained by the proposed method, which consists of three series of PFR reactors with a bypass stream.These reactors can be replaced by one PFR with the volume equal to the sum of the three PFRs volume (1793.31L).Residence time in this reactor is about 176 seconds.The maximum selectivity is 1.557.Fig.10shows variation of concentrations of C and D, and the selectivity S C/D along the length of the reactor.

F i g . 9 -
Optimum network obtained for the third case study F i g . 1 0 -Profile of concentrations of C and D, as well as S C/D along the reactor for case study three

N
o m e n c l a t u r e c Aij -Concentration of component A in j th stream in i th gene c Bij -Concentration of component B in j th stream in i th gene c k Aij -Concentration of component A in j th stream in i th gene in the output of k th CSTR available in-PFR reactor simulation c k Bij -Concentration of component B in j th stream in i th gene in the output of k th CSTR available in PFR reactor simulation c i0 -Concentration of component i in feed stream k i -Reaction rate constant Mu -Mutation rate Mu max -Maximum mutation rate Mu min -Minimum mutation rate OB b -Value of the best objective function OB w -Value of the worst objective function r A -Rate of reaction S C/A -Selectivity of component C to component A S C/D -Selectivity of component C to component D V C -Volume of CSTR V P -Volume of PFR V Ck -Volume of k th CSTR in PFR reactor simulation v ij -Flow rate of j th stream in i th gene v k ij -Flow rate of j th stream in i th gene in the output of k th CSTR in PFR reactor simulation x Ai -Conversion of component A in each CSTR x k Ai -Conversion of A component in k th CSTR in PFR reactor simulation x i -Mole fraction of i th component y i -Split and recycle ratios I n d i c e s i -number of gene j -number of stream k -number of CSTR in PFR reactor simulation p -product stream s -side or bypass stream in -inlet out -outlet Schweiger and Floudas 36) is the same, but the residence time and the volume of the reactor is considerably reduced.
F i g .7 -Maximum and average molar concentrations of B vs. no. of iteration for case study two F i g .8 -Second optimum network obtained for the second case study literature (i.e.