MINLP in Chemical Reaction Networks

Trabajo Fin de M´ aster MINLP in Chemical Reaction Networks Presented by: MaAsunci´on Jim´enez Cordero Supervisors: Dr. Rafael Blanquero Bravo, Universidad de Sevilla Dr. Emilio Carrizosa Priego, Universidad de Sevilla June 24, 2014 Contents 1 Introduction 3 1.1 Problem Definition . . . . . . . . . . . . . . . . . . . . . . . . 3 1.1.1 Stoichiometric Matrix and the Action Mass Law . . . . 3 1.2 Objective ............................. 7 1.3 Methodology ........................... 7 2 Fitting Concentrations 9 2.1 Weibull .............................. 9 2.2 Logistic .............................. 10 2.3 Mixtureofmodels......................... 11 2.4 Splines............................... 11 2.5 RBFGaussians .......................... 12 3 MINLP Formulation 20 3.1 Objective Function . . . . . . . . . . . . . . . . . . . . . . . . 21 3.1.1 Differential Method . . . . . . . . . . . . . . . . . . . . 21 3.1.2 Integral Method . . . . . . . . . . . . . . . . . . . . . . 22 3.2 Constraints ............................ 22 4 Continuous Relaxation and Projection over the Integers 30 5 Combinatorial Method 33 6 Examples 39 6.1 DataSet1............................. 40 6.2 DataSet2............................. 44 A Data Sets 56 A.1 DataSet1............................. 56 A.2 DataSet2............................. 58 2 Chapter 1 Introduction 1.1 Problem Definition A material balance for Nchemical species taking place in Lchemical reactions comprising a chemical reaction network can be written for certain classes of chemical reactors such as batch, fed-batch (semibatch) and continuous flow stirred tank reactor (CSTRs) ([19], [1], [10], [12], [18]). These balances are often written to reflect the pragmatic assumptions that the (homogeneous phase) reactor is operating isothermally, is well mixed, and that the overall density of the reaction mixture is not significantly changed by the occurrence of the chemical reactions within the reactor. The material balance expression for each species may be written in a very abstract form as (1.1): (1.1) dyi dt =fi(y1, . . . , yN, t, θ), i = 1, . . . , N where yiis the molar concentration of specie Eiat time t, and θis a multidimensional parameter. Equation (1.1) is a set of coupled ordinary differential equations (ODEs) that describe the dynamic behaviour of the reactive species due to chemical reactions, as represented by the Nreaction fluxes fi. 1.1.1 Stoichiometric Matrix and the Action Mass Law In this work we focus on a particular form of (1.1) studying the so called action mass law. In this model, the flux terms fiin (1.1) are directly linked to the stoichiometries of the Lreactions taking place and the kinetic rate terms of these reactions. To be more precise, we consider the N×Lstoichiometric 3 matrix Acontaining the stoichiometric coefficients for each of the Nspecies in the Lreactions: AN×L=   a1 1· · · aL 1 . . ..... . . a1 N· · · aL N    Here, al iis the stoichiometric coefficient of the chemical species Eiin the l−th reaction. By convention, al i<0 for a species that undergoes net consumption in a reaction (reactives), al i>0 for a species that undergoes net production (products), and al i= 0 for a species that is either not involved in the l−th reaction or has no net change in it. In the action mass law, reactives and products play very different roles. For this reason, it is convenient to express the stoichiometric matrix Aas the difference of two matrices, P, corresponding to the products, and R, corresponding to the reactives, AN×L=PN×L−RN×L=   λ1 1−ν1 1· · · λL 1−νL 1 . . ..... . . λ1 N−ν1 N· · · λL N−νL N   , where al i=λl i−νl i, i = 1, . . . , N, l = 1, . . . , L and, PN×L=   λ1 1· · · λL 1 . . ..... . . λ1 N· · · λL N   ;RN×L=   ν1 1· · · νL 1 . . ..... . . ν1 N· · · νL N    Example 1.1.1 Let us consider a hypothetical chemical reaction network comprising four reactive species E1, . . . , E4(N= 4) involved in two elementary reactions (L= 2) (1.2) E1+E2−→ E3+E4 E3−→ E4 4 The stoichiometric matrix Afor this network is A4×2=    −1 0 −1 0 1−1 1 1     , and the matrices Pand Rcorresponding to the products and reactives respectively are P4×2=    0 0 0 0 1 0 1 1     ;R4×2=    1 0 1 0 0 1 0 0     The component balances may now be conveniently expressed in terms of Aand the Lindividual reaction rates in matrix-vector form: (1.3) dy dt =Aδ(y, t, θ), where δis the L×1 dimensional vector of reaction rates, θis a set of parameters, and yis the N×1 dimensional vector of species concentrations y1,. . . , yNat time t. The Lreactions involved in δare, in general, nonlinear functions of the concentrations, y, and linear functions of what is named the rate coefficients, k. If elementary reactions are assumed, then the form of the Lrate terms in δ is determined uniquely by the reactives in each of the Lelementary equations. The law of mass action kinetics states (see [19]) that the rate of an elementary reaction may be assumed to be directly proportional to the collision frequency of the reactives, and hence the product of the reactive concentrations; i.e. for a general reaction lwith Nspecies y1, . . . , yNand associated stoichiometric coefficients νl 1,...νl N(reactives) and λl 1,...λl N(products), N X i=1 νl iyi kl −→ N X i=1 λl iyi, 5 the corresponding reaction rate δis given by δl=kl N Y i=1 yνl i i,for l= 1, . . . , L. For the reaction network in (1.2), we can add the values of the rate coefficients in k: (1.4) E1+E2 k1 −→ E3+E4 E3 k2 −→ E4 The rate vector δis in (1.5). (1.5) δ=k1y1y2 k2y3 Hence, the NODEs, for known initial conditions, yi(0) of the concentrations yi(t), i= 1,...N can be used to describe the temporal evolution of the species concentrations: dy dt =Aδ(y, t, k) =     −k1y1y2 −k1y1y2 k1y1y2−k2y3 k1y1y2+k2y3     Note that any given row of Adoes not uniquely correspond to an elementary reaction, nor does any Acorrespond to a unique set of elementary reactions. For instance, the following reaction network has the same Aas the network described by Equations (1.4): E1+E2 k1 −→ E3+E4 (1.6a) E1+E3 k2 −→ E1+E4 (1.6b) 6 1.2 Objective Once we know which are the most important issues in our problem, let us see the goal. The objective of this problem is the following: given a set of empirical concentrations measured in some time instants, find the stoichiometric coefficients of the matrix Aand the vector of rate coefficients kthat best fits to the empirical concentrations. In other words, we want to find the best set of parameters, A,k, such that the concentrations, y, obtained from (1.3) give the best fit with respect to the empirical concentrations. 1.3 Methodology Let us make a summary of the most important steps in the problem resolution. In chapters 2-5 a more detailed explanation is given. The first step is to build a function that fits the empirical concentrations in order to obtain also an approximation to the derivative of the empirical values. Let us recall that in this problem it is very important to know the concentrations but also the derivatives, since we have an ODE, (1.3), to be solved. For this purpose, different models are presented in Chapter 2. Once the empirical data are smoothed to a given fitting function, the following step is to formulate the fitting problem, keeping in mind all the chemical constraints and giving a mathematical formula to express them. It is important to note that, as will be seen later, we have to solve a Mixed Integer Non Linear Problem (MINLP), where each component of the stoichiometric matrix Acan only take integer values, and the components of the vector kare non negative real numbers. A mixed integer nonlinear program (MINLP), [4], is a problem of the following form:                min x,y f0(x, y) s.t. fj(x, y)≤0, j = 1,...J x∈Zn + y∈Rm + (1.7) 7 where nis the number of integer-constrained variables, mis the number of continuous variables, Jis the number of constraints, and finally, fj:Zn +×Rm +→R, j = 1, . . . , J are arbitrary functions mapping. In the next step, for a fixed stoichiometric matrix, A, we begin a twophase method in which the coefficient rate, kwill try to be determined. These phases, called the differential method and the integral method, consider the problems (1.8) and (1.9), respectively, where ediff and eint are two different error measures detailed in Chapter 3 and Kis the set of constraints that the rate reaction vector kmust obey: (1.8) (min kediff (y, t, A,k) s.t. k∈K (1.9) (min keint(y, t, A,k) s.t. k∈K Due to the difficulty of the MINLP, we propose to address first its continuous relaxation, and thus non integer values are obtained. Anyway, the continuous relaxation is multimodal, and the local search used, solved with the nonlinear package Minos, [16], is embedded in a multistart procedure. When the continuous relaxation problem is solved, now we search the integer matrix that best approximates to the continuous one, not just rounding the values (infeasible matrices can appear) but solving a new optimization problem. In this case, the solver CPLEX, [8], is used. Since it uses a branch and bound strategy, it gives the global optimum of the problem. Chapter 4 gives more information for this issue. We have said that the optimization in the differential and integral methods is performed for a fixed matrix A. However, the components of this matrix are also variables in our problem. The method used to find them is a combinatorial method based in the Variable Neighbourhood Search, VNS ([15]). More details about this can be seen in Chapter 5 Finally, Chapter 6 shows the results obtained with the methodology proposed for the data that appear in Appendix A. 8 Figure 2.4: Empirical Concentrations and Logistic fitting of the data in [5]. Figure 2.5: Empirical Concentrations and Mixture of 2 Weibull fitting of the data in [6]. 15 Figure 2.6: Empirical Concentrations and Mixture of 2 Weibull fitting of the data in [5]. Figure 2.7: Empirical Concentrations and Mixture of 2 Logistic fitting of the data in [6]. 16 Figure 2.8: Empirical Concentrations and Mixture of 2 Logistic fitting of the data in [5]. Figure 2.9: Empirical Concentrations and Spline fitting of the data in [6]. 17 Figure 2.10: Empirical Concentrations and Spline fitting of the data in [5]. Figure 2.11: Empirical Concentrations and RBF-Gaussians fitting of the data in [6]. 18 Figure 2.12: Empirical Concentrations and RBF-Gaussians fitting of the data in [5]. 19 Chapter 3 MINLP Formulation In this chapter we make a detailed description of the formulation of the problem we want to solve. The main objective of this problem is to obtain the chemical reaction model for the concentration dynamics of a set of species that best explains the observed concentrations at the observation instants. Then, taking as input the species concentrations at some observation instants, we try to identify the integer entries of the stoichiometric matrix, A, (1.1.1) and the rate constants, k. It is unlikely that a perfect match can be obtained, and therefore one needs to gauge the goodness of fit of the solution found. The empirical concentration, ˆyin, of specie Eiis known in a set of N0time instants. Such ˆyin should be close (ideally identical) to the theoretical concentration yin = yi(tn), obtained by solving the Cauchy problem expressed by the ODE (1.3) with the initial condition y(0) = y0. As goodness of fit, the following measure is proposed (3.1) N X i=1 N0 X n=1 win (ˆyin −yin)2, where coefficients win are weights reflecting the importance given by the user to errors at each species and time steps. 20 It is important to note that finding parameters Aand kby optimizing a goodness of fit measure like (3.1) is a hard challenge, since we face a mixed integer nonlinear programm (MINLP) in which some decision variables take values in integer numbers (the coefficients of A) and some take real values (the rates k). 3.1 Objective Function Let us see how to formulate the objective function of the problem depending on whether the differential or integral method is used. 3.1.1 Differential Method The idea of this method is to avoid solving numerically the equation (1.3) as part of the optimization of the goodness of fit function (3.1). First, a fitting method is used to obtain, from the empirical values ˆyin to the full range of positive numbers, ˆyi(t), t ≥0. The procedures tested are described in Chapter 2. In the differential method, the concentrations are assumed to be given, that is to say, the function ˆyi(t) obtained by the fitting methods are assumed to give an accurate fit to the concentrations functions. In other words, (3.2) ˆyi(t)≈yi(t),∀t, ∀i= 1, . . . , N As a rough approximation, we assume from (3.2) that the derivatives of the concentrations yare also accurately approximated by the derivatives of the approximations, ˆy, i.e., (3.3) d dt ˆyi(tn)≈d dtyi(tn),∀i= 1, . . . , N, ∀n= 1, . . . , N0. Then, the objective function, ediff (y, t, A,k) in (1.8), is used as a surrogate of the goodness of fit criterion (3.1), in which the law of mass action (1.3) is assumed: (3.4) ediff (y, t, A,k) = N X i=1 N0 X n=1 win d dt ˆyin − L X l=1 λl i−νl ikl N Y j=1 yνl j jn!!2 21 Hence, for a fix stoichiometric matrix, A, we try to find the rate reactions values, k, that minimize (3.4). This method is very fast in terms of computational cost but the results are not very sharp. This is due to the fact that we are assuming that the fitting of the derivative empirical concentration are the same as the derivative of the fitting of the data. 3.1.2 Integral Method The success of the differential method described in Section 3.1.1 strongly depends on the accuracy of the approximation of the derivative concentrations yby the derivative of the approximations, ˆy, (3.3). When this assumption is not true, then the differential method may not yield an accurate solution, but it can be used as an initial solution in a further (more expensive in terms of computational cost) phase, taking as objective function, eint(y, t, A,k) in (1.9), the natural goodness of fit (3.1). The only difference with the differential method is the following: each time the objective function is evaluated, the ODE (1.3), with the initial solution y(0) = y0must be numerically solved. 3.2 Constraints Let us explain the different constraints that appear in the optimization problem either with the objective function in (3.2), in the differential method, or (3.1), in the integral method. We can assume that there are Kdifferent chemical elements involved in species, named Hk, k = 1, . . . , K. Then, the formula of specie Eiis Ei= (H1)h1i. . . (HK)hKi , i = 1, . . . , N Since atoms are conserved in chemical reactions, we have N X i=1 hkiνl i= N X i=1 hkiλl i, k = 1, . . . , K, l = 1, . . . , L 22 The matrix (H)ki =hki is called the element matrix. The above equations can be written in a shorter way, (3.5) HA = 0 From the mass conservation in the l-th chemical reaction, we can easily deduce the second constraint N X i=1 Miλl i= N X i=1 Miνl i, l = 1, . . . , L It can also be expressed in a matrix form as (3.6) ATM= 0, where Mis the column vector of the molecular masses of species. The third and fourth constraints below force that, for a given reaction l= 1, . . . , L, the number of molecules of products and reactives, respectively, is bounded: (3.7) 1 ≤ N X i=1 λl i≤Cλ max, l = 1, . . . , L (3.8) 1 ≤ N X i=1 νl i≤Cν max, l = 1, . . . , L The following constraints imply that the stoichiometric coefficients of the products and reactives, respectively must be bounded: (3.9) 0 ≤λl i≤λmax, i = 1, . . . , N, l = 1, . . . , L (3.10) 0 ≤νl i≤νmax, i = 1, . . . , N, l = 1, . . . , L Finally, the constraint (3.11) kl≥0, l = 1, . . . , L, 23 forces the rate kto be non-negative, and (3.12) λl i, νl i∈Z, i = 1, . . . , N, l = 1, . . . , L forces the stoichiometric matrix to take integer values. Constraints (3.5)-(3.12) are nice, since they are linear, but they ignore relevant information, namely, •The minimum number of species, including reactives and products must be not smaller that a given value Nmin. •If a species Eiappears in a reaction as a reactive, it cannot appear in the same reaction as a product. •All reactions are simplified, in the sense that for a fixed reaction l, there exists a species isuch that λl i gcd(λl i, νl i, i = 1,...N)/∈N, or, νl i gcd(λl i, νl i, i = 1,...N)/∈N, where gcd denotes the great common divisor. •All species must appear in at least one reaction. In order to accommodate such constraint, the problem formulation must be enriched by adding new variables and constraints. Indeed, new binary variables are define: zν ilm =1,if νl i=m 0,otherwise , m = 1, . . . , νmax and, (3.13) zλ ilm =1,if λl i=m 0,otherwise , m = 1, . . . , λmax 24      min λl i,νl i,zλ ilm,zν ilm k(P−R)−Ark2 s.t. constraints in (3.20) k=kr In this problem, the variables are the integer stoichiometric coefficients of the matrices Pand R. For the resolution of this problem the global optimizer CPLEX, [8], is used. This method is going to be applied in some different experiments where for example, the weights win in the continuous relaxation of (3.20) are changed or the problem is solved with the formulation in zν ilm and zν ilm, or with zν ilm,zν ilm, λl iand νl i. As we have said before, this method is applied for determined experiments and in each experiment this methodology is applied for a fixed number of runs. Once we have all these integer matrices the optimization problem of differential method is applied, (3.20), keeping the best matrix of each experiment in terms of objective value. With these matrices, now we apply to the integral method optimization problem, (3.21). The one that gives the least value in the objective function will be taken as an initial stoichiometric matrix. A flow diagram of process explained above can be seen in Figure 4.1. 31 Initial Solution 1.Fitting concentrations exper=1 run=1 3. Projection over the integers 2.Cont. relax. in run according to exper 4. Differential Method run=run+1 run≤ Max run? Yes No 5.Int method to best matrix of diff. method in exper exper=exper+1 exper≤ Max exper? 6.End No Yes Figure 4.1: Flow diagram of the initial solution procedure. 32 Chapter 5 Combinatorial Method In Chapter 3 it is explained that, for a fixed stoichiometric matrix, A, the differential and the integral method can be applied. However, Ais also a decision variable, and thus further efforts are needed to explore the space of stoichiometric matrices. As initial matrix, the one obtained as in the process scheme explained in Chapter 4 is used. The following methodology is an adaptation to this problem of the Variable Neighbourhood Search algorithm method, VNS, [15]. VNS is a technique used in combinatorial optimization problems, such as the Travel Salesman Problem (TSP), [15]. The idea is starting from an initial solution and improve the objective value changing, what the authors call, the neighborhood in the search. Contrary to most other local search method VNS does not follow a trajectory, but explores increasingly distant neighbourhoods of the current incumbent solution, and jumps from there to a new one if and only if an improvement was made. A pseudocode of VNS ([13]) is given in Figure 5.1. In our problem a perturbation ˜ Aof the given solution Ais generated, and then the differential method, (3.20) is applied to obtain ˜ k. If ( ˜ A,˜ k) gives a better objective value, the process is repeated starting from such solution. Otherwise, a new perturbation is generated. This process continues until a stopping rule, related with number of iterations or running time is reached. After such stopping rule is fullfilled, if we want more accuracy in the results, the process can be repeated optimizing using the integral method, (3.21), instead of the differential one. 33 Figure 5.1: Pseudocode of VNS process. The key issue, as in all applications of VNS, is how to design the perturbations (neighbourhoods) of a given stoichiometric matrix. We define a perturbation of radius las follows: lreactions from Aare chosen at random and replaced by lfeasible reactions from a list previously calculated. To build this list where all the feasible reactions are, we first choose the number of nonzero species, 1, . . . , Cν max that will play as reactives. For each choice, there is a combination of 1, . . . , Cλ max products. It is important to observe that the species in each side must be different and the stoichiometric components are in the set {0, . . . , νmax}for the reactives and {0, . . . , λmax} for the products. In the cases, where the molecular weight vector, Mand/or the elemental matrix, Hare known, the feasible reactions must satisfy the constraints (3.6) and (3.5) respectively. Otherwise, a larger list is build. Let us show some examples of the VNS process in our particular cases. For the data set in A.1, where the matrix Hand the vector Mare known, the list of feasible reactions is not too large, as it can be seen in Table 5.1. In this case the values λmax =νmax = 2 and Cλ max =Cν max = 3. 34 1 -1 -1 -1 0 1 0 1 1 -1 0 1 -1 -1 1 0 -2 2 1 0 1 0 -1 -1 0 0 -1 -2 1 2 -1 1 0 0 1 2 -1 -2 1 -1 Table 5.1: Feasible reactions of the data set A.1. (5.1) Aini =      −1 0 −1 0−1−1 1−2 2 −120 110       With this matrix above, (5.1), obtained as Chapter 4 explains (more details of this matrix can be seen in Chapter 6), we can make a perturbation of radius 1. A column of Aini is chosen randomly, for instance we can assume, that it is the first one, i.e.: Aini =      −1 0 1 −1 1       Then we replace this column by one of the list in Table 5.1, for example, the last one: Aini =      1 0 −1 1 −1       35 Then, the perturbed stoichiometric matrix is the one that appears in (5.2) (5.2) Aini =      1 0 −1 0−1−1 −1−2 2 120 −110       Let us now explain the pseudocodes used to optimize the stoichiometric matrix Aand the rate coefficient vector kwith the differential or the integral method. ALGORITHM 1 1) Begin with an initial matrix Aini. 2) Optimize the variables in kwith the differential method and obtain an initial objective value, objini. 3) p:=1.Set Aopt := Aini and objopt := objini. 4) Until p=pmax , perturb Aopt and obtain A. 4.1)Perturb Aopt and obtain A. 4.2) Fixed A,optimize kwith the differential method. The objective value is obj. 4.3) If obj < objopt,then set Aopt := A, objopt := obj. 4.4) Let p=p+ 1. 5) Optimize kwith Aopt fixed and the differential method. Table 5.2: Algorithm 1 36 ALGORITHM 2 1) Begin with an initial matrix Aini. 2) Optimize the variables in kwith the differential method and obtain an initial objective value, objini. 3) p:=1.Set Aopt := Aini and objopt := objini. 4) Until p=pmax , perturb Aopt and obtain A. 4.1)Perturb Aopt and obtain A. 4.2) Fixed A,optimize kwith the differential method. The objective value is obj. 4.3) If obj < objopt,then set Aopt := A, objopt := obj. 4.4) Let p=p+ 1. 5) Optimize kwith Aopt fixed and the integral method. Table 5.3: Algorithm 2 ALGORITHM 3 1) Begin with an initial matrix Aini. 2) Optimize the variables in kwith the integral method and obtain an initial objective value, objini. 3) p:=1.Set Aopt := Aini and objopt := objini. 4) Until p=pmax , perturb Aopt and obtain A. 4.1)Perturb Aopt and obtain A. 4.2) Fixed A,optimize kwith the integral method. The objective value is obj. 4.3) If obj < objopt,then set Aopt := A, objopt := obj. 4.4) Let p=p+ 1. 5) Optimize kwith Aopt fixed and the differential method. Table 5.4: Algorithm 3 37 ALGORITHM 4 1) Begin with an initial matrix Aini. 2) Optimize the variables in kwith the integral method and obtain an initial objective value, objini. 3) p:=1.Set Aopt := Aini and objopt := objini. 4) Until p=pmax , perturb Aopt and obtain A. 4.1)Perturb Aopt and obtain A. 4.2) Fixed A,optimize kwith the integral method. The objective value is obj. 4.3) If obj < objopt,then set Aopt := A, objopt := obj. 4.4) Let p=p+ 1. 5) Optimize kwith Aopt fixed and the integral method. Table 5.5: Algorithm 4 38 Chapter 6 Examples In this chapter, results obtained by applying the procedure explained in Chapters 3 and 5 from the data in Appendix A are shown. It is important to note that in both data sets the differential method has been applied using the derivatives of the spline fittings. Observe that the column order in the stoichiometric matrix does not affect the solution, since the columns represent the number of the reaction and they can be sorted in the order we choose, keeping in mind that this order must correspond to the components of k. The maximum time allowed for the experiments in the combinatory using the differential method, first, and, then, the integral one has been 6 hours in each one, that is to say the programme has been running for 12 hours. For the experiments, the solver Scilab, [7], has been used. Furthermore, the parameters Cλ max, Cν max, λmax, νmax, Nmin that appear in (3.20) and (3.21) are taken as Cλ max =Cν max = 3 λmax =νmax = 2 Also, Nmin =3,in A.1 2,in A.2 39 6.1 Data Set 1 From the data of A.1, the initial stoichiometric matrix obtained, as explained in Chapter 4 and the kreaction rate vector obtained with this matrix, applying the differential method and integral method with their respective objective values are in Table 6.1. Initial stoichiometric matrix Aini =      −1 0 −1 0−1−1 1−2 2 −120 110       Int. method rate vector kini−diff =  0.1693065088 0.0000100000 0.0192887246   Diff. obj value 2.0493407576 ·10−12 Diff. method rate vector kini−int =  0.1068719346 0.0000100000 0.0188922986   Int. obj value 5.0654641711 ·10−8 Table 6.1: Initial stoichiometric matrix and reaction vector for the differential and integral method in data set A.1. Note, that the initial matrix has captured one of the three correct reactions. The plots of the concentrations obtained with the matrix in Table 6.1 and the values of the parameter vector kwith the differential and the integral method (see also Table 6.1) are shown in Figures 6.1 and 6.2 respectively. The next step of the procedure proposed is the combinatorial method just using the differential method (Chapter 5). The matrix Aobtained at the end of this step is (6.1) Acomb−diff =      −100 −1 1 −1 2−1 1 0−1 1 0 1 −1       40 Figure 6.6: Empirical and theoretical concentrations of the data in [6] after the combinatory using the integral method and with the differential method Plots of the concentration after the combinatory with differential method are shown in Figures 6.11 and 6.12. Note that, in spite of not having the correct parameters, the fitting is almost perfect in some species. The improvements along the process are shown in Table 6.8 and the plot in Figure 6.16. Finally, let us see in Table 6.9 the matrix obtained in the combinatory with the integral method and the rate reaction vector calculated with this matrix and the algorithms 3 and 4. Observe that the stoichiometric matrix obtained after the combinatorial with the integral method is exactly the same as the one obtained with the differential one given in Table 6.7. This is due to the fact, among others, of the absence of the molecular weights vector Mand the elemental matrix H. 47 Figure 6.7: Empirical and theoretical concentrations of the data in [6] after the combinatory using the integral method and with the integral method . Concentrations plots of the results of Algorithms 3 and 4 can be seen in Figures 6.14 and 6.15. They will be very similiar to Figures (6.11) and (6.12). Improvements along the process can be seen in Table 6.10 and Figure 6.16. 48 Figure 6.8: Improvements in the combinatory with the integral method for the data in [6] . Initial stoichiometric matrix Aini =      −1 1 −1−2 0 0 −1 1 −1−111 1000 1100       Int. method rate vector kini−diff =    0.0033319342 0.0000513809 0.0041930041 0.0032959845     Diff. obj value 5.6519473266 ·10−10 Diff. method rate vector kini−int =    0.0000100000 0.0169919665 0.0098788095 0.0010711874     Int. obj value 1.5335471601 ·10−5 Table 6.6: Initial stoichiometric matrix and reaction vector for the differential and integral method in data set A.2. 49 Figure 6.9: Empirical and theoretical concentrations of the data in [5] with the initial parameters and the differential method Figure 6.10: Empirical and theoretical concentrations of the data in [5] with the initial parameters and the integral method. 50 Stoichiometric matrix Acomb−diff =      0−2−1−2 −1 0 −1 1 −1011 1100 2−100       Alg 1. rate vector kAlg1=    0.0024306747 0.2199761353 0.0030956085 0.0039886964     Alg 1. obj value 1.4693823086 ·10−10 Alg 2. rate vector kAlg2=    0.0022631146 0.1672035153 0.0023254067 0.0037494357     Alg 2. obj value 4.2592913072 ·10−6 Table 6.7: Parameters for the Algorithms 1 and 2 of the Data Set A.2. Figure 6.11: Empirical and theoretical concentrations of the data in [5] after the Algorithm 1. 51 Figure 6.12: Empirical and theoretical concentrations of the data in [5] after the Algorithm 2. Figure 6.13: Improvements in the combinatory with the differential method for the data in [5]. 52 TIME (SECONDS) OBJ. VALUE 0 0.5651947327 ·10−9 3.062000000 0.4635914666 ·10−9 20.20300000 0.4055557529 ·10−9 129.9210000 0.3991671705 ·10−9 240.8430000 0.3923546641 ·10−9 249.8430000 0.3254454054 ·10−9 477.1090000 0.3134104704 ·10−9 479.2650000 0.2850598601 ·10−9 505.6090000 0.2719036922 ·10−9 526.5000000 0.2692205611 ·10−9 657.2960000 0.2417530342 ·10−9 660.4840000 0.1804711967 ·10−9 1433.921000 0.1804711966 ·10−9 3096.171000 0.1804711966 ·10−9 18903.46800 0.1758743592 ·10−9 19334.84300 0.1725465483 ·10−9 20681.45300 0.1719223189 ·10−9 21087.32800 0.1719223189 ·10−9 21400.87500 0.1469382309 ·10−9 Table 6.8: Improvements objective value in the combinatory with the differential method for the data in [5]. Stoichiometric matrix Acomb−int =      0−2−1−2 −1 0 −1 1 −1011 1100 2−100       Alg 3. rate vector kAlg3=    0.0024306746 0.2199761474 0.0030956084 0.0039886964     Alg 3. obj value 1.4693823086 ·10−10 Alg 4. rate vector kAlg2=    0.0021676152 0.2869917220 0.0023575863 0.0036784425     Alg 4. obj value 5.3409901223 ·10−6 Table 6.9: Parameters for the Algorithms 3 and 4 of the Data Set A.2. 53 Figure 6.14: Empirical and theoretical concentrations of the data in [5] after the Algorithm 3. Figure 6.15: Empirical and theoretical concentrations of the data in [5] after the Algorithm 4. 54 TIME (SECONDS) OBJ. VALUE 0 0.4259291307 ·10−5 Table 6.10: Improvements objective value in the combinatory with the integral method for the data in [5]. Figure 6.16: Improvements in the combinatory with the integral method for the data in [5]. 55 Appendix A Data Sets A.1 Data Set 1 The first set of data is taken from [6]. We have: •N= 5. •L= 3. •N0= 11. •Time steps: T=T1×11 =0 120 240 360 480 600 720 840 960 1080 1200  •Elemental matrix: H=H3×5=  40213 31212 62446  •Molecular weights vector: M=M1×5=102 18 60 32 74  •Initial concentrations: Y0=Y01×5=0.1 0.1 0 0.1 0  56