Learning R using a Chemical Reaction Engineering Book: Part 3

January 26, 2013
By

(This article was first published on Notes of a Dabbler » R, and kindly contributed to R-bloggers)

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

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):

G=-(x_1lnK_1+x_2lnK_2)+(1-x_1-x_2)lnP+(y_{I0}-x_1-x_2)ln(y_I)+(y_{B0}-x_1-x_2)ln(y_B)+(y_{P10}+x_1)ln(y_{p1})+(y_{P20}+x_2)ln(y_{P2})

The following constraints need to be satisfied for x_1 and x_2

\begin{aligned}  0 &\leq x_1 &\leq 0.5 \\  0 &\leq x_2 &\leq 0.5 \\  x_1&+x_2 &\leq 0.5 \\  \end{aligned}

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 Ax-b \geq 0

#  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 h(x) \geq 0. 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 x_1+x_2 > 0.5 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")})

gibbscontour

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

To leave a comment for the author, please follow the link and comment on his blog: Notes of a Dabbler » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.