# 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,
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)