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 sections A.2.2 and A.2.3 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.
\[
\begin{aligned}
PK_1y_Iy_B-y_{P1}=0, \
PK_2y_Iy_B-y_{P2}=0
\end{aligned}
\]
The relation between the variables \( y_I,y_B,y_{P1},y_{P2}\) and extent of reactions \( x_1,x_2\) are:
\[
\begin{aligned}
y_I=\frac{y_{I0}-x_1-x_2}{1-x_1-x_2}, \
y_B=\frac{y_{B0}-x_1-x_2}{1-x_1-x_2}, \
y_{P1}=\frac{y_{p10}+x_1}{1-x_1-x_2}, \
y_{P2}=\frac{y_{p20}+x_2}{1-x_1-x_2}
\end{aligned}
\]
Here I have used R package rootSolve for solving the above set of equations. 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.1334 0.3507
##
## $f.root
## F1 F2
## 6.162e-15 1.621e-14
##
## $iter
## [1] 7
##
## $estim.precis
## [1] 1.119e-14
##
# solution to the equations
xans$root
## [1] 0.1334 0.3507
The solution to the equations is accessed from the variable xans$root which in this case is (0.1334,0.3507)
Function Minimization
Here the function to be minimized is
\[
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\)
\[
0 \leq x_1 \leq 0.5, 0 \leq x_2 \leq 0.5, x_1+x_2 \leq 0.5
\]
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 = -3720
dg2 = -4490
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.1331 0.3509
##
## $value
## [1] -2.559
##
## $counts
## function gradient
## 100 NA
##
## $convergence
## [1] 0
##
## $message
## NULL
##
## $outer.iterations
## [1] 3
##
## $barrier.value
## [1] 8.695e-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.314
## Min(hin): 0.01602
## par: 0.1333 0.3507
## fval: -2.559
## Min(hin): 0.01598
## par: 0.1332 0.3508
## fval: -2.559
## Min(hin): 0.01598
## par: 0.1332 0.3508
## fval: -2.559
xans3
## $par
## [1] 0.1332 0.3508
##
## $value
## [1] -2.559
##
## $counts
## function gradient
## 97 14
##
## $convergence
## [1] 0
##
## $message
## NULL
##
## $outer.iterations
## [1] 4
##
## $barrier.value
## [1] 0.0008698
##
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 evaluted in the outer
# function
fcont = function(x, y) eval_f0(c(x, y))
fcontv = Vectorize(fcont, SIMPLIFY = FALSE)
z = outer(x1, x2, fcontv)
# 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 solving nonlinear equations (fsolve and 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.
This post was created with R 2.15.1 and knitr 0.7.1 (converted markdown to html)