Identification problem in plug-flow chemical reactors using the adjoint method

Identification problem in plug-flow chemical reactors using the adjoint method A. Berm´udeza, N. Estebanb, J.L. Ferr´ına,∗, J.F. Rodr´ıguez-Caloc, M.R. Sillero-Denamield aDepartamento de Matem´atica Aplicada. Universidad de Santiago de Compostela. 15782 Santiago de Compostela. Spain bInstituto Tecnol´ogico de Matem´atica Industrial. 15782 Santiago de Compostela. Spain cCentro de Tecnolog´ıa Repsol. 28935 M´ostoles (Madrid). Spain dIMUS. Universidad de Sevilla. 41012 Sevilla. Spain Abstract The aim of this work is to solve identification problems in plug-flow chemical reactors. For this purpose an adjoint-based algorithm for parameter identification problems in systems of partial differential equations is presented. The adjoint method allows us to calculate the gradient of the objective function and the constraint functions with respect to the unknown parameters significantly reducing the computer time. This leads to solve a minimization problem, in which an objective function is defined in order to quantify the mismatch between the observed data and the numerical solution of the parameterized chemical model. For solving the initial and boundaryvalue problem we use finite-difference schemes. More precisely, we propose a second-order BDF method initialized with a first-order one. The algorithm proposed was implemented in a computer program and some numerical results are shown. The efficiency of the adjoint method, compared with the classical formula of incremental quotients, is also presented. Keywords: chemical kinetics, plug-flow reactor, identification, adjoint method, integral method 1. Introduction Important models from different fields of engineering and applied sciences depend on unknown parameters which must be identified to define completely the system. Examples are common in biological systems, chemical reaction mechanisms, fluid mechanics, etc. Usually, the parameters are estimated by solving a minimization problem, where the objective function depends on the solution of the parameterized model. These models can be systems of numerical, ordinary (ODE) or partial differential equations (PDE). Many optimization algorithms consist of gradient-based iterative methods, which require, at each step, the computation of the derivatives of the objective function and possibly of the constraint functions with respect to the optimization variables. Parameter identification problems present two classical difficulties: first the number of variables and parameters of the model can be large and second the governing equations usually involve nonlinear functions of variables and parameters so, in practice, they have to be solved numerically. Therefore, a discrete minimization problem is considered instead of the continuous one, where the objective function depends on the numerical solution of the parameterized model. ∗Corresponding author Email addresses: [email protected] (A. Berm´udez), [email protected] (N. Esteban), [email protected] (J.L. Ferr´ın), [email protected] (J.F. Rodr´ıguez-Calo), [email protected] (M.R. Sillero-Denamiel) Preprint submitted to Computers & Chemical Engineering October 28, 2016 In this paper, we consider a quite general PDE-constrained optimization problem. The main goal is to propose an algorithm to compute the global solution of the discrete optimization problem obtained from the previous continuous one by solving the model of the system using some discretization procedure. More specifically, we concentrate our efforts in introducing an efficient strategy for computing the derivatives of the discrete objective function with respect to the optimization variables. Basically, three different alternatives could be used: classical finite-difference methods, automatic differentiation or adjoint method. The first one is very sensitive to the discretization step, leading to round-off errors when it is small and to truncation error when it is large. Additionally, the needed computer time is usually very high, mainly for problems with a large number of optimization variables. The second one is a set of techniques to transform a program that calculates numerical values of a function, into a program that computes the values of the derivatives of this function with about the same accuracy and efficiency as the function values themselves. Therefore, its cost for our problem will be of the same order or, very likely, higher than solving the state equation. Some references about automatic differentiation are Bartholomew-Biggs et al. (2000) and Neidinger (2010). Furthermore, there are a lot of commercial and free packages as, for example, ADOL-C, CasADi or FABDAB++. The third one comes from the mathematical theory of optimal control of systems governed by partial differential equations, as developed by J.L. Lions (see Lions (1971) and Lions (1972)). In optimal control theory, the gradient of the objective function is calculated indirectly by solving the so-called adjoint equation. Since the adjoint equation is linear, and thus of reduced complexity, the cost of obtaining the derivatives of the objective function with respect to each optimization variable from the solution of the adjoint equation is low. This procedure, which has been successfully used in different fields and applied to different problems (see, for example, Daescu et al. (2003) or Sandu et al. (2003) where it was applied in chemical kinetic systems, Jameson (1988) or Pironneau (1974) in design optimization in fluid dynamics, Maute et al. (2003) in aeroelastic systems, Meric (1985) in thermoelasticity or Vautard et al. (2000) in atmospheric chemistry), will be the one chosen in our work. In this paper, we consider models of reaction systems where the physico-chemical magnitudes depend on the particular position in the reactor (distributed systems), but only on its axial coordinate, i.e., they are supposed to be constant in the cross-sections. They correspond to the so-called plug-flow reactors (PFR). We assume that stoichiometric matrix, empirical measurements for variables and functional expressions for the reaction velocities are known, but the latter depend on a set of unknown parameters. The proposed method to estimate the values of these unknown parameters consists in minimizing the deviation between the experimental measurements and the numerical solution of the parameterized model. This optimization problem is solved by iterative methods which require, at each step, the computation of both the objective function and its gradient. In order to obtain the value of the objective function, a finite difference approximation is used to solve the model which is an initial-boundary value problem. On the other hand, we employ the adjoint method to obtain the gradient of the objective function with respect to the unknown parameters. We notice that with this method these derivatives are exactly and efficiently calculated. In addition, since nonconvex nonlinear functions appear in the optimization problem, the standard local search tools may get easily trapped in local optima. Therefore, Global Optimization techniques must be used to obtain the global optimum of the problem. In this work, the well-known Variable Neighborhood Search (VNS) method is applied. For more details on the efficiency of this methodology, the reader is referred to Carrizosa et al. (2012), Mladenovi´c and Hansen (1997) or Mladenovi´c et al. (2008). The identification problem presented in this paper is less treated in the bibliography (see, for example, Rodrigues et al. (2015a)) than the case where a stirred tank reactor is considered. In the latter, the mathematical model describing the temperature and species concentrations consists 2 of ordinary differential equations. Different methods can be found in the bibliography identifying the kinetics of the reaction system from experimental data. On one side we have the so-called simultaneous approach where all reaction rates are identified simultaneously (see, for example, Bardow and Marquardt (2004) or Bhatt et al. (2012)). On the other side we have the so-called incremental approach where the identification problem can be decoupled dealing with each reaction in an individual manner. This is mainly done using the concept of extents as, for example, in Bhatt et al. (2011), Brendel et al. (2006), Rodrigues et al. (2015b) or Srinivasan et al. (2012). Let us emphasize that the adjoint method we consider in this paper can be applied to any of the two approaches. The present paper is organized as follows. A general PDE-constrained optimization problem arising from the parameter identification problem of a reaction system in a PFR is presented in Section 2. For the sake of completeness, the model is deduced for the unsteady and non-isothermal case so the temperature is also calculated by including the energy conservation equation. The numerical methods for solving this optimization problem are introduced in Section 3. They include a BDF second-order finite-difference scheme for solving the system of partial differential equations, a so-called simultaneous integral method for reaching the global optimum and the discrete adjoint method for computing the gradient of the objective function with respect to the optimization variables. In Section 4, we apply the proposed algorithm to a plug-flow reactor and show the advantage of using the adjoint method. We finish with some conclusions. 2. Identification problem In this section we develop the identification problem for isothermal or non-isothermal PFRs when the kinetics of the reactions follow the generalized functional form: δl(θ, y,z) = Ble−El/Rθ Ml Y k=1 "N X r=1 Gl kryr+ N+Nc X r=N+1 Gl krzr−N+bkl#αkl , l = 1, . . . , L (1) where θis the temperature, yare the species concentrations vector, zare the catalysts concentrations vector, Ris the universal gas constant, Nis the number of species, Ncis the number of catalysts, Lis the number of reactions and Bl(Arrhenius pre-factors), El(activation energies), Gl kr,bkl and αkl are the unknown parameters to be found. In what follows, all these parameters to be identified will be denoted by vector Θ. We also introduce notation xfor the vector containing species concentrations and temperature, i.e., x= (yt, θ)t∈RN+1. Kinetics given by equation (1) were provided by the engineers of Repsol (an integrated oil and gas Company), taking into account that it can represent all the situations observed in the experimental setups they are trying to simulate. Then, being provided experimental data for species concentrations and temperature at the reactor exit, ( d xe,n), at several times, (Se), and for different experiments (E), the identification problem is formulated as min Θb J(Θ),(2) being b J(Θ) = J(x(Θ),Θ) := X e∈E N+1 X i=1 X te n∈Se ωien(xe,n i(Θ)−d xe,n i)2,(3) where xe,n i(Θ) is the i-th component of the solution of the mathematical model describing the PFR at time te n∈Seand at the reactor exit. Methodology to solve the optimization problem (2)-(3) is the goal of this paper. This will be addressed in the next section. 3 2.1. PFR reactor model Plug flow reactors, also called continuous tubular reactors or piston flow reactors, are reactors with cylindrical geometry. In principle, due to viscosity, the velocity of the flow in a pipe is null on the wall. Moreover, in the laminar regime, the velocity profile is parabolic with the maximum at the central axis of the pipe. However, the plug flow is a simple model where the velocity is assumed to be constant on any cross-section of the pipe (but it may depend on time). Furthermore, it will be assumed that the physico-chemical magnitudes may vary along the reactor but they are constant in any cross-section. This means that they depend on time and on one space variable. Let zbe the axial coordinate of the reactor of length L. Then z∈[0,L]. By assuming that (i) all thermodynamic magnitudes depend only on zand t, (ii) the diffusive, conductive and viscous dissipation terms can be neglected, (iii) there is no external volumetric heat source, and (iv) the density of the mixture in the reactor is constant and equal to the density of the input mixture, the full model for a non-adiabatic PFR is ∂y ∂t +v∂y ∂z =Aδ(θ, y,z),(4) ∂θ ∂t +v∂θ ∂z =1 ˆ w0(θ)·y−∆ˆ H(θ)·δ(θ, y,z) + 2g R(θext −θ),(5) y(0, t) and θ(0, t) are given,∀t∈[0, T],(6) y(z, 0) = y0(z), θ(z, 0) = θ0(z),(7) being v(t) the axial velocity of the flow at time t,Athe stoichiometric matrix, δthe reaction velocity vector with components given by (1), θext the outside temperature, gthe convective heat transfer coefficient and Rthe radius of the reactor. The components of vector ˆ w(θ)∈RNare defined by ˆwi(θ) = Miˆei(θ),(8) where Miis the molecular mass of species i, while the components of the L-dimensional vector ∆ˆ H(θ) := Atˆ w(θ) (9) are the heat of reactions at temperature θ. Moreover, the function of temperature ˆei(θ) is the specific (i.e., per unit mass) internal energy of the i-th species. It is given by ˆei(θ) := e∗ i+Zθ θ∗ ˆci(s)ds, (10) where e∗ iis the internal energy of formation of the i-th species at temperature θ∗and ci= ˆci(θ) is the specific heat of the i-th species at temperature θ. It must be taken into account that v(t) is supposed to be given (in fact, it is the velocity of the input current at the inlet of the reactor which can be obtained as its volumetric flow rate, uin(t), divided by the area of the reactor cross-section). For the sake of simplicity, we will write system (4)-(7) in the following more compact form: ∂x ∂t +v∂x ∂z =F(t, x,z),(11) x(0, t) are given,(12) x(z, 0) = y0(z) θ0(z)(13) 4 where F(t, x,z) = Aδ(θ, y,z) 1 ˆ w0(θ)·y−∆ˆ H(θ)·δ(θ, y,z) + 2g R(θext −θ)!.(14) 3. Numerical methodology 3.1. Solver As we have said at the beginning of Section 2, the initial-boundary value problem (11)-(14) must be solved to obtain xe,n(Θ) = (ye,n(Θ)t, θe,n(Θ))t. This can be done by using, for instance, a finite difference method both in time (t) and in space (z). For the sake of stability, the partial derivatives with respect to zwill be approximated by using an upwind differentiation formula, whereas they are taken implicit in time. For time discretization, we have chosen a first-order backward differentiation formula (BDF1) for the first time step and a second-order backward differentiation formula (BDF2) for the subsequent ones. Let us denote by ∆tand hthe steps for the time and spatial discretizations, respectively, and Kand Jthe number of time and spatial mesh points, named ti=i∆t,i= 0, . . . , K and zj=jh, j= 1, . . . , J. Thus, ∆t=T/K and h=L/J. By denoting xi jand zi japproximations of x(ti, zj) and z(ti, zj), respectively, and vi=v(ti), the resulting discretized equations, for j= 1, . . . , J, can be written as (for the sake of simplicity superscript ewill not be included): •For t=t1: x1 j−x0 j ∆t+v1 x1 j−x1 j−1 h=F(t1,x1 j,z1 j).(15) •For n= 2, . . . , K: 3xn+1 j−4xn j+xn−1 j 2∆t+vn+1 xn+1 j−xn+1 j−1 h=F(tn+1,xn+1 j,zn+1 j).(16) Let us notice that each time step amounts to solve Jnonlinear systems of N+1 equations (one system for each spatial mesh-point except for the leftmost one) and N+ 1 unknowns (Nspecies concentrations and temperature). They can be written in a more compact form as follows: e Bxh=e F(17) where e B=       B10 C1B2 C2C1B3 ......... 0C2C1BK        (18) is a K×Kblock-matrix and e F= (b1t,...,bKt)t.(19) 5 Bi,i= 1, . . . , K,C1and C2are J×Jblock matrices, each of them of size (N+ 1) ×(N+ 1), and bi∈RJ×(N+1). Their expressions are B1=           1 + ∆t hv1I O −∆t hv1I1 + ∆t hv1I −∆t hv1I1 + ∆t hv1I ...... O−∆t hv1I1 + ∆t hv1I           ,(20) Br=           1 + 2∆t 3hvrI O −2∆t 3hvrI1 + 2∆t 3hvrI −2∆t 3hvrI1 + 2∆t 3hvrI ...... O−2∆t 3hvrI1 + 2∆t 3hvrI           , r = 2, . . . , K, (21) C1=   −4 3I O ... O−4 3I   ,(22) C2=   1 3I O ... O1 3I   ,(23) b1=       x0 1+∆t hv1x1 0+ ∆tF(t1,x1 1,z1 1) x0 2+ ∆tF(t1,x1 2,z1 2) . . . x0 J+ ∆tF(t1,x1 J,z1 J)        ,(24) b2=       1 3x0 1+2∆t 3hv2x2 0+2 3∆tF(t2,x2 1,z2 1) 1 3x0 2+2 3∆tF(t2,x2 2,z2 2) . . . 1 3x0 J+2 3∆tF(t2,x2 J,z2 J)        ,(25) br=       2∆t 3hvrxr 0+2 3∆tF(tr,xr 1,zr 1) 2 3∆tF(tr,xr 2,zr 2) . . . 2 3∆tF(tr,xr J,zr J)        , r = 3, . . . , K, (26) where I=I(N+1)×(N+1). 6 3.2. Integral method In this subsection we are going to explain how the optimization problem (2)-(3) is solved. It must be kept in mind that there are real variables in Θand integer variables too, according to the user’s choice. Due to this, we will make a process in several steps: 1. As initial solution, the best of several solutions chosen at random is selected, thus providing a first value for the objective function. The values of the involved parameters are within user-specified ranges. (a) If all the parameters are real numbers, we start with the improvement phase, in which the continuous problem is optimized. The solution is further improved by using the VNS. (b) If there are integer parameters, a multi-start is planned for these parameters, whereas optimization is only done for the real ones. Again, the solution is further improved by using the VNS. 2. For the vector of unknowns, coming from the second step, we have to make a perturbation of radius rof each component (starting with r= 1), taking into account whether the components must take real or integer values. (a) Firstly, we consider the case where the chosen unknown takes a real value. We define an interval where the perturbation will take place. In the neighborhood of radius r, this interval has the following form: c−c−Lc rmax ∗r, c +Uc−c rmax ∗r,(27) where cis the best value found up to now for the parameter that we are considering, Lcand Ucare the lower and the upper bounds of c, respectively, and rmax (maximum value that rcan take) is the maximum number of allowed neighborhoods. In general, this interval is asymmetric. Once this interval is defined, a random number is selected in it. (b) Secondly, let us see the case where the parameter can only take an integer value. We define the same previous interval for the integer unknown, and later, we modify this interval as follows: dc−c−Lc rmax ∗re,bc+Uc−c rmax ∗rc.(28) An integer number γbetween dc−c−Lc rmax ∗reand bc+Uc−c rmax ∗rcis chosen. Here, b.c and d.edenote the floor (it returns the integer part of the number) and the ceiling (it returns the integer part of the number plus one) functions respectively. Then we proceed with a multi-start with random assignment of parameters within the intervals previously explained. If an improvement in the value of the objective function is achieved, value of ris reset to 1. Otherwise r=r+ 1. This process is repeated until either the maximum number of neighborhoods rmax or the computation time provided by the user to complete this step is reached. 3. The vector of parameters obtained in the previous step is introduced as a new initial solution for optimization, thereby providing a value of the objective function to be compared with that obtained in step 1. The process is repeated until the maximum run-time is reached or until the value of the objective function is less than a user-defined tolerance. In this optimization process, the derivatives of the objective function are needed, and for this we are going to use the adjoint method which will be explained below. 7 3.3. Adjoint method In order to obtain the gradient of the objective function involved in the optimization problem (2)-(3), we could proceed by an approximation of the partial derivatives through numerical derivation as follows ∂b J ∂Θi (Θ)≈b J(Θ+εiei)−b J(Θ) εi ,(29) for i∈ {1, . . . , Np}, where Npis the number of parameters to identify, εiis the discretization step for the i-th variable and eiis the i-th vector of the canonical basis of RNp. However, it is well-known that the computing time for this methodology is very high, due to the fact that the initial-boundary value problem (4)-(7) must be solved once for each partial derivative. Moreover, the choice of steps εiis not an easy task. We notice that, if these steps are not properly chosen, then the accuracy is very poor leading to a bad convergence of the optimization algorithm. As a less expensive and more accurate (in fact, exact) alternative, we will use the adjoint method described below. The first step is to write the continuous optimization problem (2)-(3) in a discrete form: min Θb Jd(Θ),(30) where b Jd(Θ) = Jh,∆t(xh(Θ),Θ) := X e∈E N+1 X i=1 X te n∈Se ωien(xe,n∗(n) d,i (Θ)−d xe,n i)2,(31) xhis the solution of (17) and xe,n∗(n) d,i is the i-th component of xhat time t∗ n∗(n)=te n. For the sake of simplicity, we assume that the set of observation instants is a subset of the time mesh, i.e., Se⊂ {t∗ n=n∆t, n = 0, . . . , K}. Differentiation of the discrete objective function, b J∆t, with respect to Θleads to ∇Θb Jd(Θ) = ∇Θxh· ∇xJh,∆t(xh(Θ),Θ) + ∇ΘJh,∆t(xh(Θ),Θ).(32) The gradient of the discrete variables xhwith respect to the optimization variables, ∇Θxh, could be obtained from the differentiation of the solution of system (15)-(16) by using the implicit function theorem. However, it is well-known that the computational cost of this method is high when the number of optimization variables becomes large, because one should solve Nplinearized state equations, one for each independent variable. For this reason, we are going to use the adjoint approach, that we recall below. First, by using (17) we notice that ∇Θxhsatisfies ∇Θxh·e Bt=∇Θe F.(33) Now, we introduce the adjoint state pas the solution of the following linear system: e Btp=∇xJh,∆t(xh(Θ),Θ).(34) Then, by using (33) and (34) in (32), we obtain ∇Θb Jd(Θ) in terms of the adjoint state, namely ∇Θb Jd(Θ) = ∇Θe Fp +∇ΘJh,∆t(xh,Θ).(35) Notice that this method only needs to solve one linear system, the so-called adjoint-state equation (33). The computation of the derivatives appearing in the above equations, namely, ∇xJ∆t, ∇ΘJ∆tand ∇Θe Fcan be done analytically. Remark 3.1. Because of the triangular block structure of matrix e B, the solution of the linear system in (34) can be obtained in a recursive manner, starting from the final time, by solving J×Klinear systems of size N+ 1. Moreover, the derivative of the discrete objective function can be obtained cumulatively from the adjoint state. In practice, these two steps are done simultaneously. 8 4. Numerical experiments In this section we will show the results of some numerical experiments obtained with the developed code. First, we build a synthetic experiment where the expressions of the kinetics of the reactions are provided, which allow us to solve the model of the reactor (equations (4)-(7)) with different initial values to generate the ”experimental” data. Second, we use the integral method introduced in the previous section to solve the raised identification problem and show the obtained results, comparing them with the exact ones. Finally, we remark on the difference in computing time required for calculating the gradient of the objective function using the conventional method of incremental quotients and the adjoint one. 4.1. Describing the test problem The performance of the integral method for solving the identification problem will be shown by constructing a test case consisting of a set of six reactions (L=6) involving twelve species (N=12) and one catalyst (Nc=1), and using certain given expressions for the kinetics of these reactions. More precisely, the reactions, catalyzed by C1, are E1+E2−→ E3+E4,(36) E2−→ E3+1 2E5,(37) E2−→ E6+E7,(38) E2−→ E8+E9,(39) E2−→ E10 +E11,(40) E4+E7−→ E12,(41) whereas expressions of the components of the reaction velocity vector are δ1(θ, y)=5·1010e−62737/Rθy1y2z1,(42) δ2(θ, y)=2·1010e−71097/Rθy2 2z0.5 1,(43) δ3(θ, y) = 1010e−66914/Rθy2 2z1,(44) δ4(θ, y)=5·109e−85734/Rθy2z1.1 1,(45) δ5(θ, y)=4·109e−133829/Rθy2z0.5 1,(46) δ6(θ, y)=2·109e−63987/Rθy4y7z1.2 1.(47) The reactor considered in this experimental setup has one meter long (L= 1) and a volume of 0.7854 m3, and the outside temperature is θext = 300 K. Specific heat of all the Nspecies is equal to 2000 J/kg K, whereas their molecular weights are M= (0.05,0.14,0.12,0.07,0.04,0.122,0.018,0.07,0.07,0.09,0.05,0.088)tkg/kmol.(48) Heat of the reactions at the reference temperature of 298.15 K are ∆ ˆ H(θ) = (1.5,1,0.8,0.7,0.6,0.5)t J/kmol. Finally, the overall heat transfer coefficient, g, is equal to 6 W/m2K (see (5)). With these data, 10 synthetic experiments (all of them with a catalyst concentration of 0.001 mol/l at any time) were built considering different initial and boundary values of species concentrations and temperature as shown in Table 1. Each data set corresponds to 11 equally spaced times (K= 11), until reaching the end time T= 20 seconds. The values given in this way by the mathematical model governing the reactor will be used as experimental data in the identification process, without adding perturbations. When a row in Table 1 has two lines, it means that the 9 Sandu, A., Daescu, D.N., Carmichael, G.R., 2003. Direct and adjoint sensitivity analysis of chemical kinetic systems with KPP: Part I - Theory and software tools. Atmos. Environ. 37, 5083–5096. doi:10.1016/j.atmosenv.2003.08.019. Srinivasan, S., Billeter, J., Bonvin, D., 2012. Extent-based incremental identification of reaction systems using concentration and calorimetric measurements. Chem. Eng. J. 207-208, 785–793. doi:10.1016/j.cej.2012.07.063. Vautard, R., Beekmann, M., Menut, L., 2000. Applications of adjoint modelling in atmospheric chemistry: Sensitivity and inverse modelling. Environ. Model. Softw. 15, 703–709. doi:10.1016/S1364-8152(00)00058-X. 16