Scientific paper Development of a Mathematical Model for the Dynamic Optimization of Batch Reactors, and MINLP Synthesis of Plug-flow Reactors in Complex Networks Marcel Ropotar,ab* Zdravko Kravanjab a Tanin Sevnica kemicna industrija d.d., Hermanova cesta 1, 8290 Sevnica, Slovenia b Faculty of Chemistry and Chemical Engineering, University of Maribor, P.O. Box 219, Maribor, Slovenia * Corresponding author: E-mail: marcel.ropotar@uni-mb.si Received: 28-08-2007 Dedicated to the memory of professor Vojko Ozim Abstract This paper describes the development of a robust and efficient reactor model suitable for representing batch and plug-flow reactors (PFRs) in different applications. These would range from the nonlinear (NLP) dynamic optimization of a stand-alone batch reactor up to the mixed-integer nonlinear (MINLP) synthesis of a complex reactor network in overall process schemes. Different schemes for the Orthogonal Collocation on Finite Element (OCFE) and various model formulations, in the case of MINLP model, were studied in order to increase the robustness and efficiency of the model. A deterministic model for known kinetics was obtained for batch and PFR reactors and extended for uncertainties in process parameters and reaction kinetics when the kinetics is unknown. Different variations of the developed model were applied to certain problems, as examples. The first motivating example was the dynamic optimization of batch reactor and the second the MINLP synthesis of a process scheme for the production of allyl chloride. The NLP version of the model with moving finite elements was found to be the most efficient for representing a batch reactor in the dynamic optimization example, and PFR trains in the process synthesis example. Keywords: Batch reactor, PFR reactor, orthogonal collocation, NLP, MINLP, process synthesis 1. Introduction Kinetics in batch and PFR reactors is described using differential equations with time as independent variable in the case of batch reactors and retention time, reactor length or volume in the case of PFRs. These equations represent complex optimization problems, even in small and simple examples, because equation-oriented solvers cannot handle differential equations. The use of OCFE in optimization models of batch or PFR reactors has become a well-established numerical method. It is used to convert and approximate differential equations into a set of nonlinear algebraic equations in a variety of applications, ranging from dynamic optimization of a single stand-alone batch reactor up to MINLP synthesis of complex reactor networks in overall process schemes. Over the last decade modelling, dynamic optimization, and on-line optimization have been the main research categories regarding the optimization of batch reactors. The modeling category is usually oriented towards a more realistic description of a batch reactor1 and towards the use of special modeling techniques and strategies in cases of imperfect knowledge of kinetic studies involved, e.g. the use of tendency models2 or a sequential experiment design strategy based on reinforcement lear-ning.3 The second category is related to more advanced aspects of dynamic optimization in respect of batch reactors, e.g. robust optimization of models, characterization by parametric uncertainty,4 or stochastic optimization of multimodal batch reactors.5 Finally in work relating to online optimization, which is currently the prevailing activity, different control schemes have been proposed, e.g. feedforward/state feedback laws in the presence of disturbances and nonlinear state feedback laws for batch processes with multiple manipulated inputs have been developed.6'7 Due to the complexity involved, dynamic optimization problems are regarded as difficult research tasks. On the other hand, the synthesis of reactor networks in overall process schemes is even more complex because we are dealing with discrete (selection of units, connectivity, etc.) and continuous (temperatures, flows, pressures, etc.) decisions simultaneously, which give rise to complex high-combinatorial mixed-integer nonlinear problems. Several methods have been developed for solving MINLP problems and one of the more efficient is the outer approximation (OA) algorithm8 and its extensions. It is also possible to solve MINLP reactor network synthesis problems using the geometrical approach,9 based on the attainable region (AR) theory or even by more efficient hybrid approaches which combine both methods.10 The geometrical approach, based on the AR theory, was first used for constructing an attainable region in the concentration space for 2-dimensional problems,11 and then for multi-D problems.12 Recently a novel concept of time-dependent Economic Regions (ERs) was incorporated into the MINLP synthesis of reactor networks within the overall process scheme.13 ER is obtained when economic criteria (e.g. annual profit, cost) are plotted vs. volume, residence time, or some other variable in contrast to the Concentration Attainable Region (CAR), which is constructed using technological criteria (e.g. conversion, selectivity, yield). A very important objective when optimizing reactor systems is to obtain reliable and feasible solutions, even in the presence of uncertain parameters. A lot of work has been carried out so far in design under uncertainty. For example, a novel approach was developed for the evaluation of design feasibility/flexibility, based on the principles of the deterministic global optimization algorithm a-BB14 and a two-stage algorithm for design under uncertainty and variability was proposed.15 Efficiency when solving the above-mentioned reactor optimization problems depends significantly on the method applied to solve the embedded differential-algebraic systems of equation. From among different variations of OCFE methods, the one with fixed finite elements is the most straightforward and easiest for modeling batch and PFR reactors. However, when using fixed finite elements directly it is impossible to explicitly model the optimal retention times of the batch reactors nor the optimal outlet concentrations and conditions. Consequently, the use of flexible finite elements or moving finite elements is regarded as a conventional approach for overcoming these difficulties16. This model, however, seems to be more nonlinear because the length of the final element is converted into a variable. The aim of this paper is to present the development of mathematical models suitable for optimization of batch and PFR reactors, which may either stand alone or be combined in complex reactor networks embedded within overall process flowsheets. Different schemes and strategies are applied to modelling and solving these dynamic and synthesis problems. The objective is to identify the most robust and efficient solution procedure. 2. Experimental - Numerical Procedure The following four-step procedure was proposed for solving optimization problems that contain differential-algebraic systems of equation: Simulation: During the first, optional step, simulation was performed using the MATHCAD professional package. The simulation is useful for preliminary analysis of a given kinetic system's behaviour, and to provide a good initial point for NLP or MINLP. Model formulation: During the second step, a differential-algebraic optimization problem (DAOP) model was converted into an NLP or MINLP model. Differential equations were approximated into a set of nonlinear algebraic equations by the use of OCFE, and an integral term in the objective function was approximated by the Gaussian integration formula. Solution: During the next step, either NLP or MINLP optimization was performed for the developed model. Analysis: During the last step, sensitivity analysis was carried out by one-parametric NLP or MINLP optimization with production rate (demand) as a varying (uncertain) parameter. Sensitivity analysis can be upgraded for flexible dynamic optimization where uncertain parameters are included directly in the optimization. When process synthesis is carried out, ER can be constructed during this step, with reactor volume as varying parameter. 2. 1. Dynamic Optimization of Batch Reactor Motivating example: NLP and MINLP models for the optimization of batch and PFR reactors were developed, based on a motivating example of a batch reactor (Fig. 1), where consecutive reaction A ^ B ^ C is carried out and B is the desired product. Since the reaction is endothermic, the system can be heated and/or preheated. Whenever the optimal inlet temperature is higher than defined by the user the inlet must be preheated. The kinetics of this reaction is the following: Steam Figure 1: Batch reactor. where cA, cB and cC are concentrations of A, B, and C, respectively, k0 is a pre-exponential constant, R universal gas constant, T reaction temperature, t time, and EaA and EaB are activation energies of both consecutive reactions. The corresponding (DAOP) is given as follows: (1) AtHp ^ | df pcprní\ p e P V p cP ^««»i > 0 (DAOP) where cr, cp, rr and rp denote concentration and reaction rate for reactants and products, respectively, Z profit, Nb number of batches, C cost coefficients, topt optimal reaction time, T0 inlet temperature, Tb desired temperature, ArH reaction enthalpy, cp specific heat capacity, p density of reactive mixture, #heat/cool and #preheat/precooi heat-flow for heating/cooling and heat-flow for preheating/precoo-ling, respectively. The objective is to maximize revenue for a certain number of batches. Costs for reactants and utilities are subtracted from the profit of the product sale (eq. (1)). Note that the cost function of the utility is integrated into the objective function over the whole time period. Eqs. (2) and (3) represent differential equations for the production rates of reactants and products, respectively, while eq. (4) is the differential equation for heat-flow. Heat-flow for preheating or precooling is calculated using eq. (5). The (DAOP) model for motivating example, shown above, cannot be used in equation-oriented solvers because they cannot handle differential equations. Therefore, the differential equations have to be converted into a set of nonlinear algebraic equations. In this way a (DAOP) model is converted into an NLP or MINLP model suitable for optimization. The OCFE method was applied to approximate the differential equations. Below we present three variations of the OCFE method: i) one with fixed finite elements, ii) one with moving finite elements and iii) one with fixed but partly employed finite elements. Deterministic NLP and MINLP models for dynamic optimization of batch reactors were developed, based on these variations. In order to handle deviations of uncertain parameters, the deterministic models were upgraded with flexibility constraints. Thus, the flexible dynamic optimization models were finally developed. 2. 1. 1. NLP Model Formulation Let us first describe the case of using the OCFE method with fixed and partly-employed finite elements. The following deterministic model (DFIX-NLP) was obtained, which is usually non-flexible or is flexible only for very small deviations of uncertain parameters: (6) Residual equations for component balances: h \ ^ n t /=0 k=0,k*j ' = f>Ci,;- f] - t « - U t=0.txj 1 jl tn-t, -k ■ e R T" -r + Ic ■ e R T" • r - 0 Icq ■ e ~Ea.B RT" c =0 v/=1,2,..j: }■ V/= 1,2,,..Aï (p® ~ } = iCBM -CB)+ ( V/ = \,2,...NE (15) All final elements are set as equal and total time is defined as a sum of the lengths of all finite elements: Aa, = AaM V/ = ...NE ^ V/ = 1,2,...Ms (16) (17) y,V/= (21) = 1,2,...NE (22) Figure 3: The graphical representation of 5 moving finite elements with 3 collocation points. where yl denotes binary variable for finite element l. Ineq. (19) is applied to ensure that all finite elements up to the last selected one are, in fact, selected. If the corresponding finite element is rejected, ineq. (20) forces t;opt to zero. When the element is not the last one, ineqs. (21) and (22) are applied to force the tlopt of each finite element into the upper bound. Hence, all the selected finite elements are fully exploited for integration, except the last one where the optimal time is continuously defined by the Legendre polynomial between the element's bounds. Note that, in contrast to the NLP model where integration is distributed equally and continuously within all the finite elements, here integration is only applied to the selected finite elements. In the case of NLP optimization, the number of fi- nite elements has to be set in advance and is, thus, usually overestimated in order to satisfy a given error tolerance, whereas, in MINLP cases, it is explicitly modelled in order to adjust it simultaneously to the minimal number of elements, during the optimization process. In the case of the MINLP model, the robustness of the model was studied with respect to the use of different model formulations motivated by recently developed alternative convex-hull model formulation (ACH)17 Namely a comparison was made between, Big-M formulation, conventional convex hull (CCH) and alternative convex hull formulation (ACH). In addition, the following representations of OAs in the solution point xk for the Outer Approximation/Equality Relaxation (OA/ER) algorithm were compared: Big-M formulation: h(xk)+ Vrft(x*)' a: - V Xh(xk)' xk < M(I - y) CCH representation: vAx-U^Ax^x'-hix'iv ACH representation: )' x < Vxh{xk )T xf + (v )' (x* -xf)- tfx' ))v Unlike CCH representation, where the continuous variables x are usually forced into zero values when the corresponding disjunctives are false, in ACH the variables are forced into arbitrarily-forced values, xf. Finally, in order to obtain better approximation of the OCFE method, additional inequality constraints for approximation error were included in the NLP and MINLP models (ineqs. (23)-(25)). These inequalities minimize the difference between values from the current finite element and the starting point of the next finite element 2. 1. 3. Flexible Dynamic Optimization Different changes in operating conditions, costs, quality of raw material etc. could significantly affect the steady-state operation of the process and, hence, the desired amount and quality of the product. Such changing parameters are called uncertain parameters and processes which can tolerate these changes are regarded as flexible processes. For this reason, it is important to consider uncertainty, and hence flexibility, as additional constraints when obtaining flexible process solutions. The main task of flexible design is to obtain optimally over-sized design variables for process equipment, which assure feasible solutions over the entire range of uncertain parameters using optimal investment costs. To ensure flexibility, besides nominal conditions, optimization has to be performed simultaneously at critical points, which is achieved by setting uncertain parameters at vertex points when the problem is convex. Thus, optimization at the critical vertex points serves as a flexibility constraint. In this way the deterministic model was extended by the flexibility constraints, defined at all vertex points. This was done for all equations, inequalities, and state and control variables, except for design variables because they must correspond to all vertex points simultaneously. Consequently, the size of the process equipment is valid for every possible combination of uncertain parameters. Objective function was approximated at the nominal point. The model obtained is, therefore, (NC + 1) times bigger than the deterministic model, where NC is the number of vertex points. In the case with three uncertain parameters and eight vertex points, the model is nine-times larger vs. the Gaussian integration method with five quadrature points for continuous distributions for every uncertain parameter where the model would be 133 times bigger than the deterministic Ml k=0,krj ' ft *U fl J-a k=().k* K K n t„ -1 - A if ' ki L _ m,M RTl=0J+[ (23) (23) Jm0 k=0,k*jtft hi fil fk! _L R-Ti=0J+l —Kt) ■ C (24) t £ r ^TMj+l g>g.i=Q,MI P-Cp V-P-Cp where e is an error tolerance, e.g. 10-3. (25) one. Therefore it is very useful, when no probability functions are known for uncertain parameters, to approximate objective function at a nominal point and, hence, very large model and expensive calculations are avoided. The following flexible model with moving finite elements (FMOV-NLP) was obtained where initial concentration of A, temperature, demand, pre-exponential factor and activation energy were defined as uncertain parameters: Equal distribution of finite element lengths: (31) (32) V« = \,2,...NC +1 (33) (FMOV-NLP) (26) Residual equations and component balances: ^R.ïVu (Ain ) X.cn.iiu ' ] J j=o K k=0,k*j tjh, 'tin -k0 -e Rç„u (f ilu ) - T"! CC. t!u ' J-0 hhl hh, k-QJc+J t ji„ hhi -kfy ■ e ft RTu,, ' CA,:/:/ + ' e - Bjlv RTn„ ■ cn.,„, = 0 V/=I,2,..JS: }VI = 1,2,...NE V« = 1,2,.JVC+1 Additional component balance: v" A,» — Cajiu ) — (^n.i/d — ' ri.i, )+ {ccjiii ~ ('CM ) (27) Residual equations and energy balances: (28) where additional index u is defined for NC vertex points and the nominal point. Optimal outlet point by Legendre polynomials: *=o ,k*j tju hu K K 11 _ * K K 1 _ t „opt . n 7-opt vt . n —_^ B> ¿ Ejl. 11 _ ' Xbl Zj /'« 11 j=0 k=0,k*j 'ju rku J=0 j=0 k=a,k*j'j„