Learning R Using a Chemical Reaction Engineering Book: Part 4

The links to previous parts are listed here. (Part 1, Part 2, Part 3). In this part, I tried to recreate the examples in sections A.3.1 of the computational appendix in the reaction engineering book (by Rawlings and Ekerdt).

Solving a system of ordinary differential equations

This example involves reaction (Benzene pyrolysis) in a plug flow reactor. The actual reactions happening are:

The rate of each reaction is given by:

The feed to the reactor consists of 60kmol/hr of Benzene (B). The temperature of the reactor is and the pressure is .The rate constants and equilibrium constants for this example are:

# Appendix A.3.1: Solution of Differential Equations

# Benzene pyrolysis example

# Parameters
# NBf - feed benzene rate - mol/h
# R - Universal gas constant
# T - Reactor temperature K
# P - Reactor pressure atm
# k1 - rxn1 forward rate constant L/mol.h
# k2 - rxn2 forward rate constant L/mol.h
# Keq1 - rxn1 equilibrium constant
# Keq2 - rxn2 equilibrium constant
pars=c(
NBf=60e3,
R=0.08205,
T=1033,
P=1,
k1=7e5,
k2=4e5,
Keq1=0.31,
Keq2=0.48
)

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.

# Calculation of mole fraction and conversion from extent of reaction at different volumes
extoutdf=data.frame(extout)
NBf=pars["NBf"]
extoutdf$conv=(extoutdf$ext1*2+extoutdf$ext2)/NBf
extoutdf$yB=(NBf-2*extoutdf$ext1-extoutdf$ext2)/NBf
extoutdf$yD=(extoutdf$ext1-extoutdf$ext2)/NBf
extoutdf$yT=(extoutdf$ext2)/NBf
extoutdf$yH=(extoutdf$ext1+extoutdf$ext2)/NBf

Next conversion and mole fraction is plotted as a function of reaction volume

# load library ggplot2 for plotting
library(ggplot2)
# load library reshape2 for data reshaping
library(reshape2)
# plot of conversion vs volume
ggplot(extoutdf,aes(x=time,y=conv))+geom_line()+
scale_x_continuous(breaks=seq(0,1600,by=200))+xlab('Volume (L)')+ylab('Conversion')+theme_bw(20)

# plot of mole fraction vs volume
tmp=melt(extoutdf[,c("time","yB","yD","yT","yH")],id.vars=c("time"),variable.name="moleFraction")
ggplot(tmp,aes(x=time,y=value,color=moleFraction))+geom_line()+
scale_x_continuous(breaks=seq(0,1600,by=200))+xlab('Volume (L)')+ylab('moleFraction')+theme_bw(20)

sessionInfo()
R version 2.15.1 (2012-06-22)
Platform: i386-apple-darwin9.8.0/i386 (32-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] reshape2_1.2.1  ggplot2_0.9.2.1 deSolve_1.10-4

loaded via a namespace (and not attached):
[1] colorspace_1.1-1   dichromat_1.2-4    digest_0.5.2       grid_2.15.1
[5] gtable_0.1.1       labeling_0.1       MASS_7.3-18        memoise_0.1
[9] munsell_0.3        plyr_1.7.1         proto_0.3-9.2      RColorBrewer_1.0-5
[13] scales_0.2.2       stringr_0.6.1

Learning R using a Chemical Reaction Engineering Book: Part 3

In case you missed previous parts, the links to them are listed below.

• Part 1
• Part 2

In this part, I tried to recreate the examples in section A.2.3 of the computational appendix in the reaction engineering book (by Rawlings and Ekerdt).

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

# function to be minimized
eval_f0=function(x){
dg1=-3.72e3; dg2=-4.49e3; T=400; R=1.987; P=2.5
K1=exp(-dg1/(R*T)); K2=exp(-dg2/(R*T))

yI0=0.5; yB0=0.5; yP10=0; yP20=0;
d=1-x[1]-x[2]
yI=(yI0-x[1]-x[2])/d
yB=(yB0-x[1]-x[2])/d
yP1=(yP10+x[1])/d
yP2=(yP20+x[2])/d

f=-(x[1]*log(K1)+x[2]*log(K2))+(1-x[1]-x[2])*log(P)+yI*d*log(yI)+
yB*d*log(yB)+yP1*d*log(yP1)+yP2*d*log(yP2)

return(f)
}

The constraints need to be specified in the form

#  constraint
A=matrix(c(-1,-1,1,0,-1,0,0,1,0,-1),ncol=2,byrow=TRUE)
> A
[,1] [,2]
[1,]   -1   -1
[2,]    1    0
[3,]   -1    0
[4,]    0    1
[5,]    0   -1

b=c(-0.5,0,-0.5,0,-0.5)
> b
[1] -0.5  0.0 -0.5  0.0 -0.5

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.

# initial guess
xinit=c(0.2,0.2)

# minimize function subject to bounds and constraints
xans2=constrOptim(theta=xinit, f=eval_f0, grad=NULL, ui=A, ci=b,
method = "Nelder-Mead")
> xans2
$par
[1] 0.1331157 0.3509254

$value
[1] -2.559282

$counts
function gradient
100       NA

$convergence
[1] 0

$message
NULL

$outer.iterations
[1] 3

$barrier.value
[1] 8.695209e-05

The solution can be accessed from xans2$par and is (0.1331,0.3509)

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.

Definition of constraints in the format for constrOptim.nl

# load library alabama
library(alabama)
library(numDeriv)

h_ineq=function(x){
h=rep(NA,1)
h[1]=-x[1]-x[2]+0.5
h[2]=x[1]
h[3]=-x[1]+0.5
h[4]=x[2]
h[5]=-x[2]+0.5
return(h)
}

> xans3=constrOptim.nl(par=xinit,fn=eval_f0,hin=h_ineq)
Min(hin):  0.1
par:  0.2 0.2
fval:  -2.313951
Min(hin):  0.01602445
par:  0.1332729 0.3507026
fval:  -2.559282
Min(hin):  0.01597985
par:  0.1331951 0.350825
fval:  -2.559282
Min(hin):  0.01597985
par:  0.1331951 0.350825
fval:  -2.559282

> xans3
$par
[1] 0.1331951 0.3508250

$value
[1] -2.559282

$counts
function gradient
97       14

$convergence
[1] 0

$message
NULL

$outer.iterations
[1] 4

$barrier.value
[1] 0.0008697741

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.

# filled.coutour and contour are overlaid with the minimum point (0.133,0.351)
filled.contour(x1,x2,z,xlab="x1",ylab="x2",
plot.axes={axis(1); axis(2); contour(x1,x2,z,levels=c(-3,-2.5,-2,-1.5,-1,0),vfont=c("sans serif","bold"),labcex=1,lwd=2,add=TRUE);
points(0.133,0.351,pch=15,col="blue")})

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.

> sessionInfo()
R version 2.15.1 (2012-06-22)
Platform: i386-apple-darwin9.8.0/i386 (32-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] alabama_2011.9-1  numDeriv_2012.3-1 rootSolve_1.6.2

Learning R using a Chemical Reaction Engineering Book: Part 2

In case you missed part 1, you can view it here. In this part, I tried to recreate the examples in section A.2.2 of the computational appendix in the reaction engineering book by Rawlings and Ekerdt.

Solving a nonlinear system of equations

This example involves determining reaction equilibrium conditions by solving the following system of nonlinear equations.

The relation between the variables and extent of reactions are:

Here I have used R package rootSolve for solving the above set of equations to determine and . The library is loaded and the functions to be solved are defined in the R function fns.

# load library rootSolve
library(rootSolve)

# function defining F(x)=0
fns=function(x){
K1=108; K2=284; P=2.5
yI0=0.5; yB0=0.5; yP10=0; yP20=0;
d=1-x[1]-x[2]
yI=(yI0-x[1]-x[2])/d
yB=(yB0-x[1]-x[2])/d
yP1=(yP10+x[1])/d
yP2=(yP20+x[2])/d
F1=P*K1*yI*yB-yP1
F2=P*K2*yI*yB-yP2
c(F1=F1,F2=F2)
}

Next, an initial guess of (0.2,0.2) is set for the variables and the equations are solved using the function multiroot (from package rootSolve)

# initial guess for x
xinit=c(0.2,0.2)

# solve the equations
xans=multiroot(f=fns,start=xinit)

# object returned by multiroot
> xans
$root
[1] 0.1333569 0.3506793

$f.root
F1           F2
6.161738e-15 1.620926e-14

$iter
[1] 7

$estim.precis
[1] 1.11855e-14

# solution to the equations
> xans$root
[1] 0.1333569 0.3506793
>

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.

Learning R using a Chemical Reaction Engineering Book: Part 1

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.

# stoichiometric matrix
stoi=matrix(c(0,1,0,-1,-1,1,
-1,1,1,-1,0,0,
1,0,-1,0,-1,1),
ncol=6,byrow=T)
> stoi
[,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    1    0   -1   -1    1
[2,]   -1    1    1   -1    0    0
[3,]    1    0   -1    0   -1    1

# rank of the stoichiometrix matrix
rank=qr(stoi)$rank
> rank
[1] 2

# reaction rate vector
r=c(1,2,3)

Given the reaction rate r=(r1,r2)’, the rate of change of species concentration R is given by

where is the stoichiometrix matrix

# rate of change of components
R=t(stoi)%*%r
> R
[,1]
[1,]    1
[2,]    3
[3,]   -1
[4,]   -3
[5,]   -4
[6,]    4

Example A1: Estimating reaction rates

The stoichometrix matrix is input below

# stoichiometry
stoi=matrix(c(0,1,0,-1,-1,1,
-1,1,1,-1,0,0),nrow=2,byrow=T)
> stoi
[,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    1    0   -1   -1    1
[2,]   -1    1    1   -1    0    0

# number of species and number of reactions
nspec=ncol(stoi)
nr=nrow(stoi)

> nspec
[1] 6
> nr
[1] 2

# true rxn rates
r=c(1,2)
> r
[1] 1 2

# true component rates
R=t(stoi)%*%r
> R
[,1]
[1,]   -2
[2,]    3
[3,]    2
[4,]   -3
[5,]   -1
[6,]    1

Simulate 2000 measured component rates

Add random noise (normally distributed with mean 0 and standard deviation 0.05) to true species rate vector R

## simulate 2000 noise estimates
e=matrix(0.05*rnorm(2000*nspec,0,1),nrow=2000,byrow=T)
Rmeas=matrix(rep(R,2000),ncol=nspec,byrow=T)+e

The least squares estimate of reaction rate vector is

## estimate reaction rates
rest=solve(stoi%*%t(stoi),stoi%*%t(Rmeas))

I was trying different plot features in R and applying to this data of estimated rates.
I found the following function that plots scatterplot with marginal histograms

# plotting scatterplot with histogram
# downloaded from web
#  http://www.r-bloggers.com/example-8-41-scatterplot-with-marginal-histograms/
#
scatterhist = function(x, y, xlab="", ylab=""){
par.default <- par(no.readonly=TRUE)
zones=matrix(c(2,0,1,3), ncol=2, byrow=TRUE)
layout(zones, widths=c(4/5,1/5), heights=c(1/5,4/5))
xhist = hist(x, plot=FALSE)
yhist = hist(y, plot=FALSE)
top = max(c(xhist$counts, yhist$counts))
par(mar=c(3,3,1,1))
plot(x,y)
par(mar=c(0,3,1,1))
barplot(xhist$counts, axes=FALSE, ylim=c(0, top), space=0)
par(mar=c(3,0,1,1))
barplot(yhist$counts, axes=FALSE, xlim=c(0, top), space=0, horiz=TRUE)
par(oma=c(3,3,0,0))
mtext(xlab, side=1, line=1, outer=TRUE, adj=0,
at=0.9 * (mean(x) - min(x))/(max(x)-min(x)))
mtext(ylab, side=2, line=1, outer=TRUE, adj=0,
at=(0.9 * (mean(y) - min(y))/(max(y) - min(y))))
par=par(par.default)

}

# scatter plot of reaction rates with marginal histograms
scatterhist(t(rest)[,1],t(rest)[,2],xlab="r_1",ylab="r_2")

There is library cars that has a command ‘scatterplot’ for plotting scatterplot with box plots and has several other options
In the plot below 50% and 90% ellipses are overlaid on the data

# scatter plot of reaction rates with marginal box plots
library(car)
scatterplot(t(rest)[,1],t(rest)[,2],reg.line=FALSE,smooth=FALSE,ellipse=TRUE)

I tried ggplot2 also for the same plot

# 2d contours/density with ggplot2 for reaction rates

# create dataframe of reaction rates
rest_df=data.frame(t(rest))
names(rest_df)=c("r1","r2")

library(ggplot2)
ggplot(data=rest_df,aes(x=r1,y=r2))+geom_point()

ggplot(data=rest_df,aes(x=r1,y=r2))+stat_density2d()

ggplot(data=rest_df,aes(x=r1,y=r2))+stat_density2d(aes(fill=..level..),geom="polygon")

My first blog post

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.