# MATLAB Plug flow reactor optimisation Problem

• MATLAB
• Will26040
In summary, the assignment is to find the optimal operating temperature and maximum product concentration of reactant B, assuming a constant temperature across the PFR length. The reaction is a series reaction: A → B → C (liquid phase). The code provided uses the fminunc function to find the optimal temperature and the ode15s function to solve the differential equations. After some adjustments based on feedback, the code seemed to be working correctly. However, a second part of the assignment required considering a nonlinear temperature profile, which the code did not account for. The code was modified, but there may still be some errors. The final part of the assignment involves an energy balance, which may require further adjustments to the code.
Will26040
TL;DR Summary
I am trying to create MATLAB code that will find the optimal operating temperature and maximum product concentration in a PFR. I have included my current attempt.
The assignment is to find the optimal operating temperature and maximum product concentration of reactant B, assuming a constant temperature across the PFR length. Please could someone help? thanks

the reaction is a series reaction: A → B → C (liquid phase)

Here is my current code which is providing T as the initial guess value:
function Main
x0 = [100]
OPTIONS = optimoptions('fmincon','Display','iter');
x = fminunc(@fun,x0,OPTIONS)
end

function obj = fun(x)
T = x(1);
Lspan = [0 2]; %reactor 2m long
y0 = [700 0 0];
[L,y] = ode15s(@(L,y) PFR1(L,y,T),Lspan,y0);
obj = -y(end,2)
end

function f = PFR1(L,y,T)
%series reactions Y is the concentrations of species and V PFR volume
Ca= y(1);
Cb= y(2);
Cc= y(3);

%Parameters
R = 8.31; %J.mol-1K-1
u = 0.25; %m.min-1

%rate constants
k1 = (1.37*10^(10))*exp(-75000/R*T); %min-1
k2 = (1.19*10^(17))*exp(-125000/R*T); %min-1

%Reaction Rates
ra = -k1*Ca;
rb = k1*Ca-k2*Cb
rc = k2*Cb;;

%Differential Equations
dCbdL = rb/u;
dCcdL = rc/u;

end

Two things:
- could you use a codeblock for Matlab code? That makes everything much easier to read
- what exactly do you want help with? Apparently this doesn't work? what is going wrong?

Will26040
yep sorry didnt know about code blocks:

Matlab:
function Main
x0 = [100]
OPTIONS = optimoptions('fmincon','Display','iter');
x = fminunc(@fun,x0,OPTIONS)
end

function obj = fun(x)
T = x(1);
Lspan = [0 2]; %reactor 2m long
y0 = [700 0 0];
[L,y] = ode15s(@(L,y) PFR1(L,y,T),Lspan,y0);
obj = -y(end,2)
end

function f = PFR1(L,y,T)
%series reactions Y is the concentrations of species and V PFR volume
Ca= y(1);
Cb= y(2);
Cc= y(3);

%Parameters
R = 8.31; %J.mol-1K-1
u = 0.25; %m.min-1

%rate constants
k1 = (1.37*10^(10))*exp(-75000/R*T); %min-1
k2 = (1.19*10^(17))*exp(-125000/R*T); %min-1

%Reaction Rates
ra = -k1*Ca;
rb = k1*Ca-k2*Cb
rc = k2*Cb;;

%Differential Equations
dCbdL = rb/u;
dCcdL = rc/u;

end

It doesn't seem to be working at the moment as it just gives me whatever I put in for x0 as the optimal temperature. The task is to find the optimal temperature (minimum) that gives the maximum Cb value. Am I going about it right? thanks

Have you tried the Matlab debugger? It helps in figuring if the in-between steps work. Now, I'm not a chemist, but these are the things that I would check:
- exp(-75000/R*T) (and the other exp as well) is always zero, after googling a bit seems like it should be exp(-75000/R/T). This explains the insensitivity on x0.
- Even after this, the value of the two k's are extremely small (9e-30 or so), which seems incorrect, or if not, then it will lead to all sorts of numerical trouble. Is this 75000 (activation energy?) correct?
- in the formula for k, the value of (1.37*10^(10)) seems very large, is this the molar concentration? Which you also specify in y0? there it is in the order of a few 100, not 10 billion.

So, that's where I would start, good luck!

Will26040 and BvU
Not a chemist either, but
• I second @Arjan82: ##\ \ e^{-E_a/(RT)}##
• T = 100 K is way too low a temperature to get the reaction going
• You have ra = -k1*Ca; so you don't want dCadL = -ra/u; but dCadL = +ra/u;

Last edited:
Will26040
thanks to you both for the help, I think I've sorted it with your advice.

I don't have MatLab, so I fumble with Octave (And with excel ) found about 357.5 K.
Did you have to make any other changes ?

Will26040
BvU said:
I don't have MatLab, so I fumble with Octave (And with excel ) found about 357.5 K.
Did you have to make any other changes ?
thanks that's what I got too. I didnt have to make any other changes

BvU
I now have a second part which I was wondering if anyone could help me please. You have to assume that the temperature profile varies nonlinearly with length (z) of the reactor:
##T(z)=a+bz+cz^2##
I have some code for it however I think it is wrong:

Matlab:
A=[-1 0 0 ]; B=[0];
Aeq=[]; Beq=[];
LB=[]; UB=[];
Nonlcon=[];
T = 350;
a = -1;
b = -1;
x0 = [T a b];
options = []
x = fmincon(@fun,x0,A,B,Aeq,Beq,LB,UB,Nonlcon,options)function obj = fun(x)
T = x(1);
a = x(2);
b = x(3);
zspan = [0 2]; % Reactor 2m long
y0 = [700 0 0]; % Initial values for species concentrations [Ca Cb Cc]
[z, y] = ode45(@(z,y) PFR2(z,y,T,a,b),zspan,y0);
obj = -y(end,2)
T1 = x(1)+(x(2).*z)+(x(3).*z.^(2));
plot(z,y(:,1),z,y(:,2),z,y(:,3))
xlabel('Length (m)');
ylabel('Concentration (mol/m3)');
legend('Ca', 'Cb', 'Cc');

end

function f = PFR2(z,y,T,a,b)
Ca = y(1);
Cb = y(2)
Cc = y(3);
T=T+(a*z)+(b*(z^2));

%Parameters
R = 8.31; %J.mol-1K-1
u = 0.25; %m.min-1

%rate constants
k1 = (1.37*10^(10))*exp(-75000/(R*T)); %min-1
k2 = (1.19*10^(17))*exp(-125000/(R*T)); %min-1

%Reaction Rates
ra = k1*Ca;
rb = k1*Ca-k2*Cb;
rc = k2*Cb;

%Differential Equations
%dTdz = a+(b*(z^2))
dCbdz = rb/u;
dCcdz = rc/u;

end

Im not sure if I am going about it in the right way. Should the equation for how T varies with length be a nonlinear equality constraint? thanks

Last edited:
Will26040 said:
Should the equation for how T varies with length be a nonlinear equality constraint?
Not a MATLAB (or octave) expert, but T is given and doesn't constrain anything.
But you can't have both Y and z for the same variable. I would replace T=T+(a*z)+(b*(z^2)); by T=T+(a*y)+(b*(y^2));...

##\ ##

okay thanks for all the help so far. Do you think it mostly looks right other than that?

Should be OK, yes -- but that is no guarantee

Is this going towards a gradually more detailed model with e.g. an exothermal reaction ? and then some cooling on the outside ?

##\ ##

BvU said:
Should be OK, yes -- but that is no guarantee

Is this going towards a gradually more detailed model with e.g. an exothermal reaction ? and then some cooling on the outside ?

##\ ##
Thanks. This is the last question that involves coding, however, it does say you can consider the energy balance at the end.

I think something is wrong with my code. I only seem to get one temperature in the output of x, but I am meant to be getting a temperature profile so that I can plot a graph of temperature vs length of reactor

You do the plotting inside the objective function ? That function can be called a zillion times !

##\ ##

And the temerature profile is given as far as I could see ? ##T(z)=350 - z - z^2##

Ah, complete misunderstanding at my end. What is it you now want to optimize ? Only ##z## (before it was T) it seems ?

Hi, thanks for your patience. I have attached the actual question; it still is T that I want to optimize, however you have to assume it varies non-linearly with the length- which is fixed at 2m. I have assumed the relationship
##T(z)=a+bz+cz^2##
I am not sure if I am doing it right.

Arjan82 said:
- Even after this, the value of the two k's are extremely small (9e-30 or so), which seems incorrect, or if not, then it will lead to all sorts of numerical trouble
Will26040 said:
I was just wondering in the post that you helped me with you said 'T = 100 K is way too low a temperature to get the reaction going' I was just wondering how you knew this? I think I'm meant to derive some constraints or something. thanks
I simply calculated a few k values
Inspired by previous painful experience mixing up temperatures in C and K.

In numerical solvers and such it is indeed often a good idea to restrict the range of ##T## and to work with a $$k_0' = k_0\, e^{-{Ea\over RT_0}}$$ that isn't so huge so that $$k(T) = k_0' \, e^{-{Ea\over R}({1\over T}-{1\over T_0})}$$to avoid problems with the argument of exponentiation.

##\ ##

BvU said:
Ah, complete misunderstanding at my end. What is it you now want to optimize ? Only ##z## (before it was T) it seems ?
That wasn't too good a shot either. You optimize ##x = [T, a, b] ## !

So I tried your code from #9. With a few mods it runs in Octave
Matlab:
function main %(pkg load optim can be done in command window)
A=[-1 0 0 ]; B=[0];
Aeq=[]; Beq=[];
LB=[]; UB=[];
Nonlcon=[];
T = 350;
a = -1;
b = -1;
x0 = [T a b];
options = optimset();
%x = fminunc(@fun,x0,A,B,Aeq,Beq,LB,UB,Nonlcon,options)
x = fminunc(@fun,x0,options)
end

In short: ihave no idea what A,B,Aeq,Beq,LB,UB,Nonlcon, are supposed to do.

But weven without those, it runs and terminates with a picture
and
Code:
f =

-97.503
33.805
63.698

obj = -370.26
x =

369.0975   -21.3115     6.9389

>>
which may or may not be a reasonable outcome.

I would expect you need to bring in other cosiderations to deal with question 3. After all, yield is one single number and optimizing ##a_i \ \ i \in (0...N)\ ## is N dimensional. Don't know how to tackle that.

When I look at the profile: I kind of expect ##C_b## at 2m to be at a maximum is equivalent to ##r_b =0\ ## at that point, but the figure doesn't show that.

I may find the time to plod through a book on reaction engineering at some point, but at the moment I can't help much further I'm afraid.

##\ ##

BvU said:
I simply calculated a few k values
Inspired by previous painful experience mixing up temperatures in C and K.

In numerical solvers and such it is indeed often a good idea to restrict the range of ##T## and to work with a $$k_0' = k_0\, e^{-{Ea\over RT_0}}$$ that isn't so huge so that $$k(T) = k_0' \, e^{-{Ea\over R}({1\over T}-{1\over T_0})}$$to avoid problems with the argument of exponentiation.

##\ ##That wasn't too good a shot either. You optimize ##x = [T, a, b] ## !

So I tried your code from #9. With a few mods it runs in Octave
Matlab:
function main %(pkg load optim can be done in command window)
A=[-1 0 0 ]; B=[0];
Aeq=[]; Beq=[];
LB=[]; UB=[];
Nonlcon=[];
T = 350;
a = -1;
b = -1;
x0 = [T a b];
options = optimset();
%x = fminunc(@fun,x0,A,B,Aeq,Beq,LB,UB,Nonlcon,options)
x = fminunc(@fun,x0,options)
end

In short: ihave no idea what A,B,Aeq,Beq,LB,UB,Nonlcon, are supposed to do.

But weven without those, it runs and terminates with a picture
and
Code:
f =

-97.503
33.805
63.698

obj = -370.26
x =

369.0975   -21.3115     6.9389

>>
which may or may not be a reasonable outcome.

I would expect you need to bring in other cosiderations to deal with question 3. After all, yield is one single number and optimizing ##a_i \ \ i \in (0...N)\ ## is N dimensional. Don't know how to tackle that.

When I look at the profile: I kind of expect ##C_b## at 2m to be at a maximum is equivalent to ##r_b =0\ ## at that point, but the figure doesn't show that.

I may find the time to plod through a book on reaction engineering at some point, but at the moment I can't help much further I'm afraid.

##\ ##
Thanks for all your help! that's all very useful info. Ill get back to you if I find a better solution

## 1. What is a MATLAB Plug flow reactor optimisation problem?

A MATLAB Plug flow reactor optimisation problem is a mathematical model used to determine the optimal conditions for a plug flow reactor system. This model is created using the MATLAB programming language and can help scientists and engineers design and optimize the performance of a plug flow reactor.

## 2. How is a MATLAB Plug flow reactor optimisation problem solved?

The problem is solved by using optimization algorithms and techniques, such as gradient descent, genetic algorithms, and simulated annealing. These algorithms aim to find the best values for the input parameters of the reactor system, such as flow rate, temperature, and reactor volume, that will result in the desired output.

## 3. What are the benefits of using MATLAB for plug flow reactor optimization?

MATLAB offers a user-friendly and efficient platform for solving complex mathematical problems, making it a great tool for plug flow reactor optimization. It also has a wide range of built-in functions and tools that can be used to implement and solve the optimization problem. Additionally, MATLAB has a large community of users and resources, making it easier to find help and support for your specific optimization problem.

## 4. What are the limitations of using MATLAB for plug flow reactor optimization?

One limitation of using MATLAB for plug flow reactor optimization is the need for prior knowledge and experience with the programming language. This can make it challenging for beginners to use effectively. Additionally, some optimization problems may require a large amount of computational power and memory, which can be a limitation for running the optimization model on a personal computer.

There are many online resources available for learning about MATLAB Plug flow reactor optimisation, such as tutorials, forums, and online courses. You can also refer to books and research articles on the subject. Additionally, experimenting with different optimization problems and seeking help from experienced users can also help improve your understanding of the topic.

• MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
1K
• MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
2K
• Biology and Chemistry Homework Help
Replies
1
Views
2K
• MATLAB, Maple, Mathematica, LaTeX
Replies
2
Views
10K