It looks like MATLAB, Octave and Python seem to be the preferred tools for scientific and engineering analysis (especially those involving physical models with differential equations). However as part of my learning R experience, I wanted to check out some of R tools for parameter fitting of models involving ordinary differential equations. R has packages deSolve for solving differential equations and FME for parameter fitting. The specific example here is taken from the computational appendix (A.6) of the book Chemical Reactor Analysis and Design Fundamentals by Rawlings and Ekerdt. In fact, all examples in this book are available in Octave and MATLAB.

This example involves two reactions in series

where we need to estimate the rate constants of the the two reactions and from data

The rate equations are captured in a function that is an input parameter to the ODE solver

# rate function
rxnrate=function(t,c,parms){
# rate constant passed through a list called parms
k1=parms$k1
k2=parms$k2
# c is the concentration of species
# derivatives dc/dt are computed below
r=rep(0,length(c))
r[1]=-k1*c["A"] #dcA/dt
r[2]=k1*c["A"]-k2*c["B"] #dcB/dt
r[3]=k2*c["B"] #dcC/dt
# the computed derivatives are returned as a list
# order of derivatives needs to be the same as the order of species in c
return(list(r))
}

Computing the predicted concentration for a given set of rate constants involves just solving the ODEs with intial conditions. This is illustrated for the parameter values (which is the actual parameters based on which data was generated).

This is wrapped into a function whose input is the parameters to be estimated and the output is the residuals. Here there are three concentrations that are fitted. In general, we may want to put different weights to these but in this example, they have the same weight. Also, the package FME has several nice features for fitting and uncertainty estimation. Here I have not used FME but directly tried to do the fit more to learn how this works.

Another way of estimating the uncertainty is using bootstrapping procedure. Here several simulated datasets are generated using the current estimated model and adding random normal noise to each data point (with mean=0 and variance = mean square error from model). Then the parameter is estimated for each simulated dataset. The set of parameters thus generated indicate the uncertainty in parameter estimates.

Next, the percentage of parameter estimates that fall within the ellipse is computed and found to be 93% (expected is 95%). In this case, the ellipsoidal approximation of parameter uncertainty seems adequate. But this might not be the case in general.

As I mentioned previously, the package FME has several functions for parameter fitting, MCMC simulation, sensitivity analysis but was not used in this post. Also, when running multiple simulations using ODE models, it might be better to use the C compiled version of reaction rate function. deSolve package provides that ability and describes it in this vignette. The data and R code is available here.

The governing equations for conversion versus volume in a plug flow reactor is based on extent of each of the reactions:
The initial conditions (corresponding to feed conditions ) are that the extent of reaction is zero.

The flow rates of each component along the reactor volume can be calculated from reaction extent

These are setup in a function that can be passed to an ODE solver. In this case the ODE solver we use is lsode from the R package deSolve. The inputs to the function are:

Variable over which the integration is done (Volume in this case)

The state variables of the system (Extent of the two reactions)

Parameters that are needed for description of the system (Rate constants, Temperature, Pressure, etc.)

The output from this function is the rate of change as described by the equations previously.

# load library deSolve for solving ODEs
library(deSolve)
# function that will be passed to odesolver
# vol is the variable over which the system is integrated (equivalent of time in batch reactions)
# ext is the extent of reactions 1 and 2
# params are the parameters passed to the system
rxnrate=function(vol,ext,params) {
with(as.list(c(ext,params)),{
NB=NBf-2*ext1-ext2
ND=ext1-ext2
NH=ext1+ext2
NT=ext2
Q=NBf*R*T/P
cB=NB/Q
cD=ND/Q
cT=NT/Q
cH=NH/Q
dext1=k1*(cB*cB-cD*cH/Keq1)
dext2=k2*(cB*cD-cT*cH/Keq2)
return(list(c(dext1=dext1,dext2=dext2)))
})
}

Since the reaction start only after the feed enters the reactor, the extent of reaction is zero for both reactions at the beginning of the reactor (V=0L). The set of volumes where the concentration and reaction extent is computed is chosen in this case to be from 0L to 1600L at every 50L. The ODE solver lsode from deSolve package is used to solve the system of equations.

# initial extent of reaction (zero in this case for both reactions)
extinit=c(ext1=0,ext2=0)
# Volumes where the concentration is reported (in this case 0 to 1600L at every 50L)
vols=seq(0,1600,length=50)
# Solution of the set of differential equations using lsode solver in deSolve package
extout=lsode(times=vols,y=extinit,func=rxnrate,parms=pars)

extout contains the extent of reaction vs volume data. That is used to compute mole fraction and conversion at different volumes along the reactor.

Function Minimization

In part 2, the reaction equilibrium conditions was determined by solving a set of nonlinear equations. Alternately, it could be determined by minimizing the Gibbs free energy of the system. The function to be minimized is (refer to part2 for information on what the variables refer to):

The following constraints need to be satisfied for and

First I used constrOptim function for minimization. We need to specify the function to be minimized

Next, the function is minimized using constrOptim (starting from an initial guess of (0.2,0.2)). Here Nelder-Mead method is used since BFGS method requires specifying the gradient of the function by the user. R taskview optimization lists other options.

Next, I also tried function minimization with ConstrOptim.nl from the package alabama. Here the constraints are specified in terms of . Even if gradient is not supplied, this function will estimate it using finite-differences.

The solution can be accessed from xans3$par and is (0.1332,0.3508).
Since this is a 2 dimensional problem, the solution can also be visualized using a contour plot of the function.

# Region of interest: 0.01<=x1<=0.49, 0.01<=x2<=0.49
x1=seq(0.01,0.49,by=0.01)
x2=seq(0.01,0.49,by=0.01)
# vectorizing function eval_f0 so that it can be evaluated in the outer function
fcont=function(x,y) eval_f0(c(x,y))
fcontv=Vectorize(fcont,SIMPLIFY=FALSE)
> z=outer(x1,x2,fcontv)
There were 50 or more warnings (use warnings() to see the first 50)

The warnings are produced in regions where since the function is not defined in those regions. This is not an issue for the contour plot (it will just ignore those regions) we will plot next.

MATLAB/Octave functions for function minimization (fmincon) have been used in Chemical Engineering computations for a long time and are robust. R has traditionally not been used in this domain. So it is hard to say how the functions I have used in this blog will perform across the range of problems encountered in Reaction Engineering.

The solution to the equations is accessed from the variable xans$root which in this case is (0.1334,0.3507)

MATLAB/Octave functions for solving nonlinear equations (fsolve) have been used in Chemical Engineering computations for a long time and are robust. R has traditionally not been used in this domain. So it is hard to say how the functions I have used in this blog will perform across the range of problems encountered in Reaction Engineering.

Chemical Reactor Analysis and Design Fundamentals by J.B. Rawlings and J. G. Ekerdt is a textbook for studying Chemical Reaction Engineering. The popular open source package Octave has its origins to the reaction engineering course offered by Prof. Rawlings. This book is accompanied by Octave and Matlab code for solving typical problems encountered in Reaction Engineering.

I figured that maybe one way to learn R is so see whether I can code some of th examples from this book in R. I am by no means suggesting that R can replace MATLAB/Octave for engineering problems but merely it is a way for me to learn the language.

I started with the computational appendix listed in the book’s website and am trying to work through some of the examples there. It will be good to refer to the computational appendix to follow the R code below.

I am interested in math, computations and visualizations. I am also interested in learning R and Python which are popular open source and free programming languages. There are big communities around these languages. As I learn from others I wanted to have a place to catalog my learnings and so have created this blog for that purpose.