Two Optimization Methods to Determine the Rate Constants of a Complex Chemical Reaction Using FORTRAN and MATLAB

: Problem statement: For chemical reactions, the determination of the rate constants is both very difficult and a time consuming process. The aim of this research was to develop computer programs for determining the rate constants for the general form of any complex reaction at a certain temperature. The development of such program can be very helpful in the control of industrial processes as well as in the study of the reaction mechanisms. Determination of the accurate values of the rate constants would help in establishing the optimum conditions of reactor design including pressure, temperature and other parameters of the chemical reaction. Approach: From the experimental concentration-time data, initial values of rate constants were calculated. Experimental data encountered several types of errors, including temperature variation, impurities in the reactants and human errors. Simulations of a second order consecutive irreversible chemical reaction of the saponification of diethyl ester were presented as an example of the complex reactions. The rate equations (system of simultaneous differential equations) of the reaction were solved to get the analytical concentration versus time profiles. The simulation results were compared with experimental results at each measured point. All deviations between experimental and calculated values were squared and summed up to form a new function. This function was fed into a minimizer routine that gave the optimal rate constants. Two optimization techniques were developed using FORTRAN and MATLAB for accurately determining the rate constants of the reaction at certain temperature from the experimental data. Results: Results showed that the two proposed programs were very efficient, fast and accurate tools to determine the true rate constants of the reaction with less 1% error. The use of the MATLAB embedded subroutines for simultaneously solving the differential equations and minimization of the error function was very fast in solving such problems, as compared to the FORTRAN program, which, although resulting in fast and accurate results, yet, requiring the use of a library of external subroutines. Conclusion: Any of the two proposed methodologies could be used to determine the rate constants of any complex reaction at a certain temperature. The proposed programs were independent of the nature of the reaction, only the rate equations and the initial conditions had to be modified for any new reaction.


INTRODUCTION
Chemical kinetics: When different chemical compounds are present under appropriate conditions, the chemical reaction will occur. Kinetics is concerned with the dynamics of chemical reactions such as the way by which reactions take place and the rate (velocity) of the process. The keystone for a reaction mechanism is its rate law. It describes the relation between the velocity of a reaction and the concentration of chemical reactants (Green and Perry, 2007).

Rate constant:
The rate constant is defined as the rate of concentration of a substance involved in the reaction with a minus or plus sign attached, depending on whether the substance is a reactant or a product. The rate constant is the constant of proportionality "k" in the rate equation as expressed by Eq. 1: Where: rate = mol dm −3 s −1 k = The rate constant C A and C B = Concentrations in mol dm −3 a and b = Orders of reaction with respect to A and B respectively The rate constant "k" is defined using Arrhenius equation (Smith, 1983) and is given by Eq. 2: Where: A = The frequency factor E A = the activation energy R = The gas constant T = The absolute temperature in Kelvin Determination of the value of the rate constant is of great importance, since it helps to determine the value of the rate of any reaction applying its rate equation. In the industrial applications of kinetics, knowledge of chemical rate equation is essential in establishing the optimum conditions of pressure, temperature, feed composition, space velocity, extent of conversion and recycling (Smith, 1983).
Surprisingly, the rate constant is not actually a true constant! It varies, for example, if we change the temperature of the reaction, add a catalyst, or change the catalyst. The rate constant is constant for a given reaction only if all we are changing is the concentration of the reactants (Santen and Niemantsverdriet, 1995).
The units of the rate constant "k" depend on the overall order of reaction ((a+b) in Eq. 1) and is given by Eq. 3 (Smith, 1983): Consecutive reactions: The reactions in which the reactants are converted to products through one or more intermediate stages are called consecutive reactions (Upadhyay, 2006). The overall reaction is a result of several consecutive steps. Every stage has its own reactants and rate constant. A consecutive reaction of the second order is given by Eq. 4a and b: Numeric solution of the rate equation: The dynamics of the elementary reactions can mathematically be expressed by a set of coupled differential equations that needs to be integrated simultaneously in order to derive the reactants' concentrations for the entire reaction. Unfortunately, only very basic rate laws have explicit mathematical solutions, i.e., solutions that can be written in algebraic form. In general, there is no analytical solution of the corresponding system of ordinary differential equations. These have to be solved numerically. Many methods are available for performing numerical integration of differential equations. Experience with a wide range of typical chemical engineering problems had shown that, with almost no exceptions, the fourth-order numerical integration method (Runge-Kutta method) is very convenient for kinetic problems. This method has the advantages of being self starting, not requiring integration at equally spaced values of time and being readily available as a library routine at most computer packages (Stoer and Bulirsch, 2002).
Simulation tools: Several programming languages and computer-aided mathematical solvers had appeared with the goal of solving complex mathematical problems (Castillo et al., 2001). FORTRAN is one such programming language. It has the advantage that several library subroutines and programs are readily available for use (Dorn and McCracken, 1972;Franks, 1972).
Modern trends for determination of the rate constants depend on new simulation tools instead of laboratory experimental data. Further optimization techniques are implemented to minimize the error between theoretical (computed) results and experimental data.
In (Arx et al., 1998), a genetic algorithm technique is used for finding a global minimum for the error function. In (Aloisio and Francisco, 2000), a simulation is done for a reaction that requires high binding energy and the rate constants are calculated from theoretical data affected by noise. In (Bijlsma and Smilde, 2000), instead of using the popular curve fitting technique to extract rate constants from the reaction equations, curve resolution methods adapted to estimate reaction rate constants from spectral data are used. In (Bijlsma et al., 2001), the same method is implemented using Ultra Violet visible (UV-vis) spectroscopy. In (Tongcheng et al., 2005), Numeric Genetic Algorithm (NGA) and Tabu Search (TS) are united to optimize the estimated rate constants of a consecutive reaction. In (Zhao et al., 2006), a nonlinear least squares regression is used to fit the kinetic profiles. Genetic algorithm is adopted to optimize the parameters of the postulated mechanisms. In (Kłos et al., 2008), the rate constants of the removal of hydroxides OH are theoretically calculated for a vibrationally excited molecule and atomic oxygen. In (Shafiei et al., 2008), the rate constant of the electrochemical oxidation of acetaminophen is estimated by comparing the experimental cyclic voltammograms with the digital simulated results using Monte Carlo algorithm. In (Ponton and Bartlett, 2009), the United States' Air Forces conducted a research work for studying the thermal rate constants for gas phase reactions. Such quantities are critical input to atmospheric and combustion models. They also provide information to designing procedures to destroy stock piles of nerve agents and other waste products. It is often very difficult to experimentally measure such rate constants, making 'predictive' quantum chemistry often the best route to their determination. Unlike experiment, rate constants can also be computed under severe conditions including very high temperatures.

Optimization techniques:
A typical engineering problem can be represented as follows: A process can be represented by some equations or perhaps solely by experimental data. The goal of optimization is to find the values of the variables in the process that yield the best value of the performance criterion (minimum cost, maximum profit, minimum losses, maximum use of resources, minimum effort and minimum error). The describing factors of the process (equations or experimental data) and the performance criterion constitute the optimization problem (Edgar et al., 2001;Venketaraman, 2009).
Optimization is the third stage in the process of designing a system. The first stage is modeling, or simulating, the system. The second step is to decide what is to be optimized, that is, to construct the socalled objective function. This function may be of single variable or of multi-variables, may be linear or non-linear and finally may be constrained by certain conditions or unconstrained. Optimization implies either maximizing or minimizing the objective function (Bak et al., 2002).
From elementary calculus, the function has a horizontal tangent at the point where a maximum or a minimum occurs. For a function of one variable this means that the derivative is zero at that point, whereas for a function of many variables it means that all the partial derivatives at that point are zero. In both cases, equating the derivative to zero is not sufficient; it is only necessary and we may find some other global maximum or minimum in the neighborhood. Moreover, determining the type of the point (maximum or minimum) requires higher derivatives (Corsano et al., 2009).

Aim of the work:
The aim of this research work is to develop two computer programs for determining the rate constants for the general form of any consecutive second order irreversible reaction at a certain temperature, as an example of a complex reaction, using FORTRAN and MATLAB. The development of such programs can be very helpful in the control of industrial processes as well as in the study of the reaction mechanisms.

MATERIALS AND METHODS
Part I: Getting the initial rate constants from the experimental data: The saponification of di-esters is taken as a case study to verify the validity of the proposed programs. This reaction is a consecutive second order irreversible reaction at a certain temperature. The kinetic equations for the reaction are given by Eq. 5a and b: Or in a symbolic form: These rate equations are modeled and the resultant set of simultaneous differential equations is given in Eq. 7a-e: • The first step in determining the rate constants of the reaction is to get the concentration-time data from laboratory experiment • The experimental concentration-time data are plotted using either FORTRAN or MATLAB (or any similar program, like Excel) • For any arbitrary time instant, get the slopes of the plotted curves • Substitute the slopes and the concentrations at the chosen time instant in the rate Eq. 7a-e, to get the rate constants of the reaction k 1 and k 2 , which are the approximated or initial values of the rate constants Part II: Optimizing the rate constants: Using FORTRAN: Complications of computing higher order partial derivatives to determine the minimum of the objective function had resulted in derivative-free optimization methods. Rosenbrock (1960); Hooke and Jeeves (1961) and Powell (1964) have introduced several algorithms for finding the minimum of multivariable, unconstrained, nonlinear functions with the advantage that no derivatives are required. These algorithms are found readily as FORTRAN subroutines. In this study, A FORTRAN program developed by Rosenbrock (1960) is used to determine the rate constants of the reaction that result in minimum error (objective) function (Kuester and Mize, 1985). The algorithm proceeds as follows: 1. Simulation of the kinetic reactions is performed with the aid of a FORTRAN program. The simulation program reads the initial values of the rate constants, the initial concentration of the reactants and the simulation start, end and interval times. The summation starts from the initial time and ends at the final time. The value of the objective function is stored.
3. The first variable k 1 is stepped a distance s 1 and the new value of the function F is evaluated. 4. If the value of F decreased, the move is termed a success and s 1 is increased by a factor α, where α≥1.
If the value of F increased, the move is termed a failure and s1 is decreased by a factor β, where 0< β≤1 and the direction of movement is reversed. 5. The second variable, k 2 , is in turn stepped a distance s 2 . The same acceleration or deceleration and reversal procedure is followed in consecutive repetitive sequences until a success (decrease in F) or failure (increase in F) have been encountered. 6. The procedure terminates when the convergence criteria is satisfied, as shown in Fig. 1.
Using MATLAB: 1. Simulation of the kinetic reactions is performed with the aid of a MATLAB computer program (mfile). The simulation program reads the initial values of the rate constants, the initial concentration of the reactants and the simulation start, end and interval times. The summation starts from the initial time and ends at the final time. The value of the objective function is stored. 3. Then we vary k 1 and k 2 in all directions. This results in eight alternatives: (i) Increase k 1 and increase k 2 (ii) Increase k 1 and keep k 2 constant (iii) Increase k 1 and decrease k 2 (iv) Keep k 1 constant and increase k 2 (v) Keep k 1 constant and decrease k 2 (vi) Decrease k 1 and increase k 2 (vii) Decrease k 1 and keep k 2 constant (viii) Decrease k 1 and decrease k 2 For each case, recalculate the concentration-time data and eight new values for F are calculated. The rate constants corresponding to the minimum F are stored and considered improved rate constants for final or next iteration. 4. The iterations proceed until the absolute difference between two successive objective functions is less than a predefined tolerance ε. The final obtained values of k 1 and k 2 will be the optimal rate constants as shown in Fig. 2.

RESULTS
Part I: Getting the initial rate constants: • Experimental data of concentration versus time are given in Table 1. Note that at time t = 0, the initial concentrations of the reactants CA and CB are:  • The concentration results are plotted versus time. C A decreases gradually whereas C B decreases more and more since it is involved in the two reactions. C C begins to increase as a result of the first reaction and decreases as it is a reactant in the second reaction. Figure 3 presents this relationship • Using MATLAB toolbox, curve fitting of the resulting plots is used to get the functions that best fit to the given experimental data. These functions are differentiated to get the slopes dCA/dt and dC B /dt. At any time instant, say at t = 10 min, get the slopes of C A and C B . These values are found to be: Table 1, the corresponding concentrations of C A , C B and C C at time t = 10 min are: Substituting these values and the slopes calculated from the previous values into the kinetic reactions equations, we get the values of the rate constants k 1 and k 2 , which are considered to be the initial or approximate (experimental) values. These values are found to be: k 1 = 5 k 2 = 2.8 Part II: Optimizing the rate constants: Using FORTRAN: The final obtained values of k 1 and k 2 , which are the optimal rate constants, are found to be: k 1 = 5.6 k 2 = 3.65 Using MATLAB: The final obtained values of k 1 and k 2 , which are the optimal rate constants, are found to be: k 1 = 5.6 k 2 = 3.64

DISCUSSION
The main differences between the two optimization methods are the number of the parameters used in calculations and the mathematical procedure of calculation, which may affect the accuracy of the results. To test the accuracy of any proposed optimization subroutine, the initial rate constant values (k 1 , k 2 ) were changed as follows: In FORTRAN, k 1 was taken as 20% bigger or smaller then the initial value, while k 2 was kept unchanged. Then the value of k 2 was increased or decreased by 20% while k 1 was kept constant. Then both k 1 and k 2 values were increased or decreased by 20%, as shown in Table 2. In MATLAB, Fig. 4 and 5 show the variation of the rate constants from the initial values to their final optimum values. Here, we go further by changing the initial rate constants by decreasing k 1 100% and increasing k 2 100%. The results are shown in Fig. 6 and 7.

CONCLUSION
In this study, two methods for determining the rate constants of a consecutive second order irreversible kinetic reaction are proposed. Both methods start by obtaining a good estimation of the rate constants by using curve fitting of the experimental concentrationtime data. These initial values, together with the initial concentrations, are fed to numerical integration programs, which utilizes the fourth order Runge-Kutta method to solve the system of simultaneous ordinary equations representing the kinetic reaction. Using least squares method, the deviations between the experimental data and the calculated data are minimized by two different methods, using FORTRAN and MATLAB to obtain the optimal values of the rate constants. The results obtained showed that both procedures are very efficient in determination of the rate constants in terms of accuracy and execution time. The MATLAB toolbox with its extended library of mathematical subroutines for numeric integration and optimization is favored to the FORTRAN method, which requires long user-built subroutines for numeric integration and optimization.
The experience of the programmer plays an important role in having accurate results. For example, the selected reactant in determining the rate constants should be involved in the rate equations containing all the rate constants to be optimized. The conversion criteria, the tolerance value and the integration step size determine to a great extent both the accuracy and the execution time of the program.
The proposed procedures could be modified to simulate any kinetic reaction, of any order and nature. Determination of accurate rate constants is very helpful in establishing the optimum conditions of pressure, temperature, feed composition, reactor design and other chemical parameters.