Format And Type Fonts CCHHEEMMIICCAALL EENNGGIINNEEEERRIINNGG TTRRAANNSSAACCTTIIOONNSS VOL. 39, 2014 A publication of The Italian Association of Chemical Engineering www.aidic.it/cet Guest Editors: Petar Sabev Varbanov, Jiří Jaromír Klemeš, Peng Yen Liew, Jun Yow Yong Copyright © 2014, AIDIC Servizi S.r.l., ISBN 978-88-95608-30-3; ISSN 2283-9216 DOI: 10.3303/CET1439118 Please cite this article as: Tóth L.R., Torgyik T., Paor D., Nagy L., 2014, Evaluation of the behaviour of objective functions in the optimization of a batch process for biodiesel production, Chemical Engineering Transactions, 39, 703-708 DOI:10.3303/CET1439118 703 Evaluation of the Behaviour of Objective Functions in the Optimization of a Batch Process for Biodiesel Production László Richárd Tóth*, Tamás Torgyik, Dávid Paor, Lajos Nagy University of Pannonia, Department of Process Engineering, 10 Egyetem utca, H-8200 Veszprém, Hungary tothl@fmt.uni-pannon.hu Modeling and optimisation of batch processes has been an emerging research field in the recent years. This work introduces the model based optimisation of a batch process, the production of fatty acid methyl esters (biodiesel). The aim of this work is to point out the importance of objective function formulation that greatly affects the optimal values of decision variables. A transesterification process takes place in a stirred vessel with jacket heating and cooling. The energy flow is manipulated by the flow rate and inlet temperature of the heat transfer liquid. A controller is also included to ensure that the temperature is near the optimal trajectory. The three studied objective functions are the maximal product concentration (purity of the product), the maximal production rate, and maximal profit generated during production. The most important variables, that may affect the process, are the initial concentrations and the temperature profile during the process. In this work mainly the latter has been studied. The use of a controller reduces the number of decision variables: instead of numerous heating/cooling energy flow values the most important points of the temperature set-point trajectory are given. They were chosen as follows: the starting high temperature, the time of switching from high to low temperature, and the time of finishing the process. There is no need for separate temperature set-point values for the heating and temperature keeping phases. The controller parameters are also included in the decision variables. All the decision variables have lower and upper limits, and there is a nonlinear constraint that represents the cooling of the reactor under a prescribed temperature at the end of the process. Numeric solvers were stuck in local minima, thus it became an interesting thought to map the objective function values throughout the whole range of the decision variables, and find a starting point closer to the optimal one. As a result we obtained approximations of the Pareto-fronts, as all the studied objective functions proved to be competing. Also the results show significant differences in the optimal temperature trajectories. The values in the case of the economical objective lie between the values in the cases of the maximal production rate and maximal concentration. 1. Introduction Process optimisation has long been an important and still emerging field of research. This work is done by modelling the process by mathematical models, formulating an objective function, gathering the constraints, and applying a numeric solver that either minimizes or maximizes the objective function by manipulating the decision variables. The optimisation task of a batch process may contain more decision variables than its continuous counterpart, because starting and end conditions highly influence the main purpose of the project. Bonvin et al. (2003) proposed an efficient optimisation method for this problem. Batch processes mainly appear in fine chemicals industry, such as in pharmaceutical, polymer or food and beverage industry, however biodiesel can also be produced in a batch reactor system. Batch processes are typically used when production volumes are low, but continuous stirred tank reactors can provide higher amount of products as well. The economic feasibility of biodiesel not only depends on the characteristic of the process (batch/continuous), but also on other factors such as different type of feedstock, transportation and storage cost of the material (Diwekar et al., 2013). 704 Researchers have made investigations with plant-based fuels, plant oils and fats to produce bio-fuel. Biodiesel, produced from vegetable oils, can be an option to replace diesel. Methyl esters of fatty acids are produced by the transesterification of vegetable oils. Triglycerides react with methanol in the presence of a catalyst at an optimal temperature range of 315-330 K (Diwekar et al., 2012). The results of Hassim et al. (2013) show that from a sustainability point of view the base catalyzed reaction is a favorable way. Yet works mainly focused either on the modeling or the numeric solver step of this procedure. In this work we would like to point out the importance of objective function selection. Formulating the objective function is often not thorough, and not explained. On the other hand it is a cliché that different objective functions have their optima at different values of decision variables. Usually it is not an easy task to implement the engineering sense into a mathematical formula. 2. Modeling of the system The purpose of modeling is to calculate the production rate of fatty acid methyl esters (FAME). This involves component balances to all components that take part in reactions in connection with FAME, and thermal modeling of the reactor, as the temperature is one of the most important factors affecting the reaction rate. We modeled a batch reactor that can be split into two main parts: the jacket and the reaction vessel. Walls, the mixer and nozzles are not contained in this model, as mainly thermal effects are studied on this level, and their contribution can be neglected. In the jacket we made some assumptions: Both heating and cooling is carried out by the same liquid (so-called monofluid). The jacket is filled totally with liquid, its volume is constant. There is no chemical reaction in the jacket. The liquid is perfectly mixed. The material properties are assumed to be constant despite the temperature change. Heat transfer towards the environment is neglected (good insulation assumed). Heat transfer through the reactor wall is described with an overall heat transfer coefficient that remains constant during the process. The resulting heat balance equation is the following (list of symbols and abbreviations can be found in Table 3.):     jr jpjj jinj j jj TT cV A TT V F dt dT  , ,   (1) The reaction vessel is treated in a similar way, but there appears the chemical reaction as heat source. We also have to construct component balance equations. We made some assumptions during the reaction vessel modeling: The reaction takes place in a single liquid phase, in the presence of an inert solvent. The gas phase does not take part in thermal energy and component transport. The liquid is perfectly mixed. The material properties are assumed to be constant despite the temperature change. The resulting heat balance equation is the following:       n i iir rpr rj rprr r rH c TT cV A dt dT 1 , ,, 1   (2) The component mass balance equations count only with the chemical reaction as a major effect. Here a matrix representation of the equation is presented: r dt dc  (3) The kinetics of the process are taken from Diwekar et al. (2012). The following reactions are taken into consideration (number over the double arrow denotes forward reaction, below the reverse reaction): FAMEDGMTG  1 2 (4) FAMEMGMDG  3 4 (5) FAMEGMMG  5 6 (6) The reaction rates are calculated by Arrhenius-type equations, and using the concentrations of the reagents. The reaction rate coefficient is calculated by Eq(1), parameters are presented in Table 1. 705 Table 1: Parameters of the chemical reactions. Cells with background colouring make up  T Reaction Component stoichiometry No. k0 (L/mol/min) Ea (J/mol) Hr (J/mol) Rate equation TG DG MG G M FAME 1 3.92e7 6,614.83 26,180 MTG cck 1 -1 1 0 0 -1 1 2 5.77e5 4,997.98 -26,180 FAMEDG cck 2 1 -1 0 0 1 -1 3 5.88e12 9,993.96 26,180 MDG cck 3 0 -1 1 0 -1 1 4 1.98e9 7,366.64 -26,180 FAMEMG cck 4 0 1 -1 0 1 -1 5 5.35e3 3,231.18 -64,580 MMG cck 5 0 0 -1 1 -1 1 6 2.15e4 4,824.87 64,580 FAMEG cck 6 0 0 1 -1 1 -1          r ia ii RT E kk , ,0 exp i=1…6 (7) As a part of our model we have to mention the controller. A PI controller has been chosen for this task, because it is sufficient in controlling a second order system like the jacketed batch reactor, while it has only two parameters. The controlled variable is the reactor temperature, the manipulated variable is a virtual input of a split-range controller. If the value of the manipulated variable is below 0.5 (50 %), then the system is in cooling mode, thus the inlet temperature of the jacket is chosen to be the lower value (293 K). Over 0.5 the heating mode is active, and the jacket inlet temperature is high (373 K). The flow rate is zero at 0.5, and changes linearly with the distance from 0.5, reaching its maximal value at 0 and 1. 3. Optimisation problem formulation There were three objectives for the process formulated: maximum product concentration, maximum rate of production, maximum profit. These functions are presented in Eq(8-10). Economical objective needs further explanation. The factors with the greatest impact on the objective value are the income from products, the cost of raw materials and energy. Separation and recirculation steps are not included in this study, and thus assumed to be ideal, although there are energy and cost efficient separation methods, like the one presented in Kiss and Ignat (2012). The needed energy is calculated by the ideal assumption that all the energy goes into the heating of the reaction mixture at start, and the heat transfer fluid during the process. Only heating is calculated, we did not go below ambient temperature, and the water, which was used for heat transfer, causes negligible costs.   endFAMEC tcJ  (8)   end endFAMEr P t tcV J  (9)        end tend ambinjjjpjambrrprrheatrawrawodprodr E t dtTTFcTTcVPcPcPV J            0 ,,,Pr 0  (10) The model of the system is used in the calculation of the objective and the constraint functions. More specifically a simulation is performed with the given design parameters, and the results are the series of the state variables (concentrations and the temperatures), and some other supplementary variables (flow rate in the jacket). The same set of simulation results can be evaluated using different objective functions. The decision variables describe the temperature set-point profile. We have the starting high temperature, the final temperature, the time of change between the two and the time for the whole process as decision variables. Also there is the starting temperature of the liquid, which is needed to model preheating prior to mixing reagents. The energy cost of preheating is included in the objective function, but the time needed for this step is not included, as it may go parallel with the reaction of the previous batch in another vessel. Finally, the PI parameters are also decision variables. We must note that the tuning of parameters is not necessarily optimal from a control perspective, but optimised according to the actual objective function. 706 All the decision variables have lower and upper limits. In some cases this is a natural limit, for example temperatures cannot be set higher than the high value of the jacket inlet temperature. In other cases it is an arbitrarily chosen value that represents some common sense consideration. If the optimal value of a decision variable had been found on an arbitrary limit, then the limit was extended to see if there exists a real optimum further, or the objective function approaches towards infinity. The problem also has nonlinear constraints, specifically the general upper limit of the reactor temperature and an upper limit to the final value of the reactor temperature. These constraints are evaluated using the simulated temperature, which is the nonlinear function of the decision variables. We needed solver methods that can handle nonlinear constraints. 4. Solution methods Three approaches were used to study the optimality and the behaviour of the objective functions. The first one is not a real optimiser algorithm; it simply evaluates the objective functions at equally distributed points between the upper and lower limit of the decision variables. As there are multiple decision variables, the points are selected to include all the possible variations of the values. This method is not computationally efficient in finding the optimal point, but rather a good way to study the behaviour of the objective functions. Also it is a good tool to select a near optimal starting point for other numeric optimiser methods. The second method for optimisation was the Matlab built-in fmincon with active-set algorithm. The results of this algorithm were quite dependent on the starting point, and had the tendency to get stuck in the local minima. By selecting a good starting point based on the previous approach, the results improved significantly. The limits were narrowed to let the optimiser search in a smaller area, thus avoiding some potentially false local minimal values. The third approach was a genetic algorithm called CMA-ES (Hansen, 2006). The results provided by this method are almost independent of the starting point. Because it is a stochastic method, it tends to get stuck less in local minima. In the case of this method some penalty functions were added to the objective functions to ensure that the optimiser respects nonlinear constraints. 5. Results and discussion The first method gave us some pictures of the objective functions evaluated at many points. To make an understandable picture, we plotted each objective functions versus the others. The resulting charts can be seen on Figure 1. One may identify that all the objective functions are competing with each other in some regions. It is also noteworthy that those cases, which do not respect nonlinear constraints, can be otherwise far better in terms of optimality. Table 2 shows the different values of the objective functions in all the three cases of optimisation. We can notice that the stochastic CMA-ES was able to reach better optimal values in all cases. In the case of optimisation to reach maximal concentration shows a serious fall in the production rate. Also this causes the economical objective to suffer. The other end is when the goal of optimisation is maximal production rate. This decreases product concentration, and depending on the solver method, the results can mean either increase or decrease in profitability. The third case is when the economic optimum is sought. In this case the other two objectives stand between their previous values. Figure 2 shows the found optimal temperature profiles. The optimisation for maximal production rate results the shortest batch. By decreasing the batch time as far as possible, the amount of product per unit time increases. The main goal of the optimiser in this case is to reach the prescribed ending temperature to satisfy the constraint. The optimisation for economic objective drives the process longer on high temperature, which generates more valuable product. On the other hand the total length of the process is still short because the Table 2: Optimal objective function values when optimisation is done by three different objective functions Value of… Algorithm Case optimised to have maximal… concentration product rate profit Concentration (mol/dm 3 ) Active-set 0.8501 0.6988 0.8060 CMA-ES 0.8535 0.5975 0.7791 Product rate (mol/s) Active-set 0.1220 0.2302 0.2126 CMA-ES 0.1152 0.2583 0.2083 Profit (HUF/h) Active-set 14,201 18,931 24,254 CMA-ES 13,124 10,050 24,361 707 faster we produce the valuable material, the faster we generate profit. The longest of all scenarios is the maximal concentration scenario. The most important factor in increasing the time is the heating in the beginning of the process. While the two other cases used the advantage of preheating, in this case preheated material would cause loss in the final concentration of the product, although it has only a weak effect. 0.2 0.4 0.6 0.8 0 0.1 0.2 0.3 0.4 0.5 0.6 Product concetration (mol/dm 3 ) P ro d u c t ra te ( m o l/ s ) -8 -6 -4 -2 0 2 x 10 4 0.1 0.2 0.3 0.4 0.5 Profit (HUF/h) P ro d u c t ra te ( m o l/ s ) 0.2 0.4 0.6 0.8 -8 -6 -4 -2 0 2 x 10 4 Concentration (mol/dm 3 ) P ro fi t (F t/ h ) Outside constraints Satisfies constraints Figure 1: Objective functions versus each other 0 500 1000 1500 2000 2500 3000 3500 290 295 300 305 310 315 320 325 330 335 340 Time (s) T e m p e ra tu re ( K ) T r (J c optimized) SP (J c optimized) T r (J P optimized) SP (J P optimized) T r (J E optimized) SP (J E optimized) Figure 2: Different optimal temperature trajectories 708 6. Conclusions It has been shown that goals of maximum production rate, maximum final product concentration and economical optimality are not satisfied in the same time. The two technological objective functions result in greatly differing optimal temperature trajectories, while the economical optimality lies somewhere between the two. This reinforces our thought that maintaining economical production needs the compromise between technological objectives. The advantage of our approach is that it uses few decision variables, and that it makes the change of the batch time possible. In further studies the problem formulation can be more refined, and by the usage of more sophisticated optimisers the increasing number of variables can be handled. Acknowledgements This work was supported by the European Union in the frame of TAMOP-4.2.2/A-11/1/KONV-2012-0071 project. Table 3: List of symbols and abbreviations A Heat transfer area 2.9 m 2 tend Ending time of batch s c Vector of concentrations mol/m 3 Tamb Ambient temperature K cp,j Specific heat capacity in jacket 4,192 J/(kg K) Tj Jacket temperature K cp,r Specific heat capacity in reactor 2,265 J/(kg K) Tj,in Jacket inlet temperature K Hr Enthalpy of reaction J/mol Tr Reactor temperature K Ea Activation energy J/mol Vj Jacket volume 0.155 m 3 Fj Flow rate of heat transfer liquid m 3 /s Vr Reactive vessel volume 0.5 m 3 JC Concentration objective function mol/m 3  Overall heat transfer coefficient 255 W/(m 2 K) JP Production rate obj. function mol/s  Stoichiometry matrix JE Economical obj. function HUF/s j Density in jacket 981 kg/m 3 k Reaction rate coefficient m 3 /(mol s) r Density in reactor 753 kg/m 3 k0 Preexponential factor m 3 /(mol s) TG Triglyceride n Number of reactions 6 DG Diglyceride Pprod Product prices HUF/mol MG Monoglyceride Praw Raw material prices HUF/mol G Gylcerol Pheat Price of energy HUF/W M Methanol r Reaction rate (vector) mol/(m 3 s) FAME Fatty Acid Methyl Ester R Universal gas constant 8.314 J/(mol K) References Bonvin D., Srinivasan B., Palanki S., 2003, Dynamic optimization of batch processes: I. Characterization of the nominal solution, Computers & Chemical Engineering, 27(1), 1-26. Diwekar U., Benavides P.T., 2012, Optimal control of biodiesel production in a batch reactor Part I: Deterministic control, Fuel, 94, 211–217. Diwekar U., Benavides P.T., Salazar J., 2013, Economic Comparison of Continuous and Batch Production of Biodiesel Using Soybean Oil, Environmental Progress & Sustainable Energy, 32(1), 11–24. Hansen N., 2006, The CMA Evolution Strategy: A Comparing Review, in J.A. Lozano, P. Larrañga, I. Inza and E. Bengoetxea (eds.). Towards a new evolutionary computation. Advances in estimation of distribution algorithms, 75-102, Springer Hassim M.H., Liew W.H., Ng D.K.S., 2013, Screening of sustainable biodiesel production pathways during process research and development (R&D) stage using fuzzy optimization, Chemical Engineering Transactions, 35, 1075-1080 DOI:10.3303/CET1335179 Kiss A.A., Ignat R.M. 2012, Energy efficient recovery of methanol and glycerol in biodiesel production, Chemical Engineering Transactions, 29, 1141-1146, DOI: 10.3303/CET1229191