Reactor Optimization with fmincon!

MIT 10.37 students and other chemical engineers in reactor enegineering might find this post interesting. Here I’m going to talk about how to apply fmincon to typical reactor engineering optimization problems.

Basic fmincon usage

Consider the optimization problem below,

it only takes three steps to translate into matlab code as shown below.

% step1: define objective function
fun = @(x)100*(x(2)-x(1)^2)^2 + (1-x(1))^2; 

% step2: set up constraints
A = [1,2];
b = 4;
Aeq = [2, 1];
beq = 1;
lb = [0,0];
ub = [0.4, 0.4];

% step3: provide initial guess and call fmincon
x0 = [1,1];
x = fmincon(fun, x0, A, b, Aeq, beq,lb, ub)

There’s an alternative way of translating constraints, instead of using A, b, Aeq, beq, lb, and ub, we can use c and ceq.

% step1: define objective function
fun = @(x)100*(x(2)-x(1)^2)^2 + (1-x(1))^2; 

% step2: set up constraints
function [c,ceq] = constraints(x)
c(1) = x(1) + 2x(2) - 4;
c(2) = x(1) - 0.4;
c(3) = x(2) - 0.4;
c(4) = -x(1);
c(5) = -x(2);
ceq = 2x(1) + x(2) - 1;
return

A = [];
b = [];
Aeq = [];
beq = [];
lb = [];
ub = [];

% step3: provide initial guess and call fmincon
x0 = [1,1];
x = fmincon(fun, x0, A, b, Aeq, beq,lb, ub, @constraints)

Matlab tutorial has well covered typical ways of using fmincon to perform optimization. Please check it out if you want to know more.

Reactor optimization example

Now with the basic knowledge of using fmincon, let’s apply it to our reactor design domain. A good example is 2016-02-22’s 10.37 recitation problem. It has a typical reaction pathway where reactant goes to desired product, which however will react to form an undesired product. As chemical engineers, we should manipulate knobs, e.g., temperature, residence time, etc., in order to maximize molar yiled of desired product.

Define objective function

As we discussed, here we want to maximize molar yield of D (diglyceride) by tuning \(T\) (temperature), \(C_{T,0}\), \(C_{CH_3O,0}\), \(C_{CH_3OH,0}\) and \(\tau\). Putting it into mathematic term would be

Translating it to matlab, we can define in below:

function obj = objFun(x)
% x(1) = tau
% x(2) = temperature
% x(3) = [T]_0
% x(4) = [CH3O]_0
% x(5) = [CH3OH]_0
% x(6) = [T]
% x(7) = [CH3O]
% x(8) = [CH3OH]
% x(9) = [D]
% x(10) = k (the kinetic coefficient)

obj = -x(9)/x(3); % molar yield of D
end

Set up constraints

Almost all the real-life chemical engineering optimization problems have certain number of constraints. They are usually from mass balance, kinetics rules, physical bounds, etc.

• Mass balance

For convenience, we assume volumetric flowrate is not changed from inlet to outlet. so we have the following equations for major species’ mass balance in CSTR.

• Kinetics rules

Kinetics coefficient k has temperature dependence according to Arrhenius law.

• Physical bounds

Physical bounds are those constraints that seem obvious to you but can largely reduce the search space for Matlab. So in order to get correct answer more quickly, we should provide as much extra information of physics as possible.

For instance, the concentrations and residence time can not be negative, and temperature should always be grater than 273 K (otherwise frozen). Here we will add another constraint where temperature cannot go beyound 350 C because of heating power limitation.

In summary we have,

Put all these constraints together in Matlab using c and ceq,

function [c, ceq] = allcon(x)
% params
% assign some easy numbers for
% Arrhenius expression params
A = exp(2);
R = 8.314;
Ea = R*600;

% ineqn
c(1) = -x(1);
c(2) = -x(2) + 273;
c(3) = -x(3);
c(4) = -x(4);
c(5) = -x(5);
c(6) = -x(6);
c(7) = -x(7);
c(8) = -x(8);
c(9) = -x(9);
c(10) = x(2) - 623;

% eqn
ceq(1) = x(10) - A*exp(-Ea/R/x(2));
ceq(2) = -x(9) + (3*x(6) - 2*x(9))*x(10)*x(7)*x(1);
ceq(3) = x(3) - x(6) - 3*x(6)*x(10)*x(7)*x(1);
ceq(4) = x(4) - x(7);
ceq(5) = x(5) - x(8) - 3*x(6)*x(10)*x(7)*x(1) - 2*x(9)*x(10)*x(7)*x(1);
end

call fmincon to optimize

Simply run the optimization by calling fmincon and feeding objFun and allcon we defined early.

% give the following initial guess for the optimal solution
% tau = 1
% temperature = 300K
% [T]0 = 4
% [C3H0]0 = 1
% [CH3OH]0 = 6
% [T] = 1
% [CH3O] = 1
% [CH3OH] = 1
% [D] = 1
% k = 1
x0 = [1,300,4,1,6,1,1,1,1,1];
[x,fval] = fmincon(@objFun, x0, [], [], [], [], [], [], @allcon)

Eventually you will have optimal yield of diglyceride to be 30%. You can find the Matlab codes here.