Optimization of Plug Flow Reactor by Varying Temperature and Pressure

I'm trying to optimize a complex reaction system to have the largest output of desired product C by varying the temperature and pressure of the reaction in an adiabatic PFR. I've created the design equation of the reaction:

Function dYdV = Designequation(v,y,T,P) 
Fa = y(1); 
Fb = y(2); 
Fc = y(3);
Fd = y(4);
% Explicit equations
T = T(1); 
P = P(1);
Cto = P / (8.314 * 10^-5) / T; 
E1 = 15000; 
E2 = 17500;
Ft = Fa + Fb + Fc + Fd; 
Ca = Cto * Fa / Ft;
Cb = Cto * Fb/Ft;
Cc = Cto * Fc/Ft;
Cd = Cto * Fd/Ft;
k1 = 0.075 * exp(E1 /1.987 * (1/300 - (1/T)));
k2 = 0.0015 * exp(E2 / 1.987 * (1 / 300 - (1 / T))); 
Fao = 50;  
ra = -k1*Ca*Cb;
rb = -k2*Cb*Cc;
% Differential equations
dFadV = ra; 
dFbdV = (2*ra)+rb; 
dFcdV = rb-ra;
dFddV = -rb;
dYdV = [dFadV; dFbdV; dFcdV; dFddV];

The design equation is solved for volume and the initial molar flow rates of A and B using ode113 with an initial temperature of 300K and 1 atm:

Vspan = [0 2]; % Range for the volume of the reactor 
y0 = [50; 25; 0; 0]; % Initial values for the dependent variables i.e Fa,Fb,and Fc, Fd
T = 300;
P = 1;
[v, y] =ode113(@(v,y) Designequation(v,y,T,P),Vspan,y0);
dFadT = y(:,1);
dFbdT =y(:,2);
dFcdT =y(:,3);
dFddT = y(:,4);
plot(v,dFcdT); hold on

Now I want to find the largest value of the molar flow rate of C (dFcdT), by varying the temperature and pressure, but I haven't found a way that I can implement that based on searches of the MatLab forum.

I apologize if this has been answered in other forum posts, but I'm desperate.

[Torsten]

So you want to maximize dFcdT(end) ? Or integral_{v=0}^{v=2} dFcdT dv ? Or something else ?

[Sean Powers]

I want to maximice dFcdT(end) while keeping the vspan 0:2 and y0 the same.

Answer 1

[Torsten]

Although I guess that there are bounds on P and T, here is a code that will give you a start.

In case you want to constrain P and T, use "fmincon" instead of "fminsearch".

function main
  format long
  T = 300;
  P = 1;
  x0 = [T P];
  options = optimset('MaxFunEvals',100000,'MaxIter',100000) 
  x = fminsearch(@fun,x0,options)
end

function obj = fun(x)
  T = x(1);
  P = x(2);
  Vspan = [0 2]; % Range for the volume of the reactor 
  y0 = [50; 25; 0; 0]; % Initial values for the dependent variables i.e Fa,Fb,and Fc, Fd
  [v, y] =ode15s(@(v,y) Designequation(v,y,T,P),Vspan,y0);
  obj = -y(end,3)
end

function dYdV = Designequation(v,y,T,P) 
  Fa = y(1); 
  Fb = y(2); 
  Fc = y(3);
  Fd = y(4);
  % Explicit equations
  Cto = P / (8.314 * 10^-5) / T; 
  E1 = 15000; 
  E2 = 17500;
  Ft = Fa + Fb + Fc + Fd; 
  Ca = Cto * Fa / Ft;
  Cb = Cto * Fb/Ft;
  Cc = Cto * Fc/Ft;
  Cd = Cto * Fd/Ft;
  k1 = 0.075 * exp(E1 /1.987 * (1/300 - (1/T)));
  k2 = 0.0015 * exp(E2 / 1.987 * (1 / 300 - (1 / T))); 
  Fao = 50;  
  ra = -k1*Ca*Cb;
  rb = -k2*Cb*Cc;
  % Differential equations
  dFadV = ra; 
  dFbdV = (2*ra)+rb; 
  dFcdV = rb-ra;
  dFddV = -rb;
  dYdV = [dFadV; dFbdV; dFcdV; dFddV];
end

[Sean Powers]

I'm extremely bad with the syntax of matlab since I've been trying to learn it here and there over the past week with little luck. Could you help me constrain the problem to temperatures between 300 and 400 and pressures between 1 and 2?

Thanks!

[Torsten]

function main
  format long
  lower_T = 300;
  upper_T = 400;
  lower_P = 1.0;
  upper_P = 2.0;
  x0 = [0 0];
  options = optimset('MaxFunEvals',100000,'MaxIter',100000,'TolFun',1e-8,'TolX',1e-8) 
  x = fminsearch(@(x)fun(x,lower_T,upper_T,lower_P,upper_P),x0,options);
  value_of_c = -fun(x,lower_T,upper_T,lower_P,upper_P)
  T = lower_T + sin(x(1))^2*(upper_T - lower_T)
  P = lower_P + sin(x(2))^2*(upper_P - lower_P)
end

function obj = fun(x,lower_T,upper_T,lower_P,upper_P)
  T = lower_T + sin(x(1))^2*(upper_T - lower_T);
  P = lower_P + sin(x(2))^2*(upper_P - lower_P);  
  Vspan = [0 2]; % Range for the volume of the reactor 
  y0 = [50; 25; 0; 0]; % Initial values for the dependent variables i.e Fa,Fb,and Fc, Fd
  options = odeset('RelTol',1e-8,'AbsTol',1e-8); 
  [v, y] = ode15s(@(v,y) Designequation(v,y,T,P),Vspan,y0,options);
  obj = -y(end,3)
end

function dYdV = Designequation(v,y,T,P) 
  Fa = y(1); 
  Fb = y(2); 
  Fc = y(3);
  Fd = y(4);
  % Explicit equations
  Cto = P / (8.314 * 10^-5) / T; 
  E1 = 15000; 
  E2 = 17500;
  Ft = Fa + Fb + Fc + Fd; 
  Ca = Cto * Fa / Ft;
  Cb = Cto * Fb/Ft;
  Cc = Cto * Fc/Ft;
  Cd = Cto * Fd/Ft;
  k1 = 0.075 * exp(E1 /1.987 * (1/300 - (1/T)));
  k2 = 0.0015 * exp(E2 / 1.987 * (1 / 300 - (1 / T))); 
  Fao = 50;  
  ra = -k1*Ca*Cb;
  rb = -k2*Cb*Cc;
  % Differential equations
  dFadV = ra; 
  dFbdV = (2*ra)+rb; 
  dFcdV = rb-ra;
  dFddV = -rb;
  dYdV = [dFadV; dFbdV; dFcdV; dFddV];
end

[Sean Powers]

Thank you for all the help already. The last part of my problem requires me to make 2-D graphs of C vs T and C vs P. I would also like to make a 3-D graph of C vs T and P.

If you had time to help with that, it would be awesome, but I understand if you are sick and tired of my needy-ness. I really wish that I had a class focused on Matlab, but unfortunately I did not and I don't have the time or skillset to learn the syntax before this problem is due.

Thanks again for all the help you've already provided.

[Bosnian Kingdom]

Did you solve the problem. Did you maximize product FC? I need to solve similar problem and i need help.

[Sean Powers]

Yes I did. I ended up evaluating the function over many temperatures and pressures using nested for loops and then finding the temperature and pressure that corresponded to the maximum value of C. Here's my code:

function iterate(Tmin,Tmax,Pmin,Pmax,Titerate,Piterate)
clc
format long
stepT = ((Tmax - Tmin)/Titerate);
stepP = ((Pmax - Pmin)/Piterate);
C = [];
T1 = []; %Kelvin
P1 = []; %bar
if(stepT == 0)
    C = zeros(1,Piterate+1);
    T1 = zeros(1,Piterate+1);  
    P1 = zeros(1,Piterate+1);
else if(stepP == 0)
    C = zeros(Titerate+1,1);
    T1 = zeros(Titerate+1,1);
    P1 = zeros(Titerate+1,1);
else
    C = zeros(Titerate+1,Piterate+1);
    T1 = zeros(Titerate+1,Piterate+1);
    P1 = zeros(Titerate+1,Piterate+1);
    end
end
TempIndex = 1;
PresIndex = 1;
if (stepT == 0)
    disp('No T iteration')
    for P = Pmin:stepP:Pmax
        T = Tmax;
        C(TempIndex,PresIndex) = optimize(T,P);
        T1(TempIndex,PresIndex) = T;
        P1(TempIndex,PresIndex) = P;
        PresIndex = PresIndex + 1;
    end
    else if (stepP == 0)
        disp('No P iteration')
         for T = Tmin:stepT:Tmax
             P = Pmax;
             C(TempIndex,PresIndex) = optimize(T,P);
             T1(TempIndex,PresIndex) = T;
             P1(TempIndex,PresIndex) = P;
             TempIndex = TempIndex + 1;
         end
        else
            disp('Regular Iteration')
            for T = Tmin:stepT:Tmax
                for P = Pmin:stepP:Pmax
                    C(TempIndex,PresIndex) = optimize(T,P);
                    T1(TempIndex,PresIndex) = T;
                    P1(TempIndex,PresIndex) = P;
                    PresIndex = PresIndex + 1;
                end
                PresIndex = 1;
                TempIndex = TempIndex + 1;
            end
        end
end
maxC = max(C(:))
minC = min(C(:))
[row, column] = find(C == maxC)
maxP = P1(row,column)
maxT = T1(row,column)
if (stepT == 0)
    plot(P1,C)
    xlabel('Pressure')
    ylabel('Output of C')
else if(stepP == 0)
    plot(T1,C)
    xlabel('Temperature')
    ylabel('Output of C')
else
    surf(T1,P1,C)
    xlabel('Temperature')
    ylabel('Pressure')
    zlabel('Output of C')
    xlim([Tmin Tmax])
    ylim([Pmin Pmax])
    zlim([minC maxC])
end
end
end

function C = optimize(T,P)
Vspan = [0 2]; % Range for the volume of the reactor 
y0 = [50/60; 25/60; 0; 0]; % Initial values for the dependent variables i.e Fa,Fb,and Fc, Fd
[v, y] =ode113(@(v,y) Designequation(v,y,T,P),Vspan,y0);
%dFadT = y(:,1);
%dFbdT =y(:,2);
C = y(end,3);
%dFcdT = y(end,3);
%dFddT = y(:,4);
end


function dYdV = Designequation(v,y,T,P) 
Fa = y(1); 
Fb = y(2); 
Fc = y(3);
Fd = y(4);
% Explicit equations
T = T(1); 
P = P(1);
Cto = P / (8.314 * 10^-5) / T; 
E1 = 15000; 
E2 = 17500;
Ft = Fa + Fb + Fc + Fd; 
Ca = Cto * Fa / Ft;
Cb = Cto * Fb/Ft;
Cc = Cto * Fc/Ft;
Cd = Cto * Fd/Ft;
k1 = 0.075 * exp(E1 /1.987 * (1/300 - (1/T)));
k2 = 0.0015 * exp(E2 / 1.987 * (1 / 300 - (1 / T))); 
Fao = 50;  
ra = -k1*Ca*Cb;
rb = -k2*Cb*Cc;
% Differential equations
dFadV = ra; 
dFbdV = (2*ra)+rb; 
dFcdV = rb-ra;
dFddV = -rb;
dYdV = [dFadV; dFbdV; dFcdV; dFddV];
end

[Bosnian Kingdom]

Thank you for answer. You write three function, but i cant't start function, where is your script?

I need to solve similar problem, but i don't have to change pressure or temperature, i only need to maximize product C from following equations:

d(FA)/d(V) = rA

d(FB)/d(V) = rB

d(FC)/d(V) = rC

d(FD)/d(V) = rD

d(FE)/d(V) = rE

d(FF)/d(V) = rF

rA = -(rA1 + rA2)

rB = -(rB1 + rB2 + rB3)

rC = rC1 + rC3

rD = rD1 + rD2 + rD3

rE = rE2 + rE3

etc.

If you have time to help me, i will write the whole code. I appreciate your help, i really need to solve this.

[Sean Powers]

Modify the function Designequation given your differential equations and reaction rates and then run iterate given a temperature and pressure range and the amount of steps you want to take in those ranges. i.e.

function dYdV = Designequation(v,y,T,P) 
Fa = y(1); 
Fb = y(2); 
Fc = y(3);
Fd = y(4);
% Explicit equations
T = T(1); 
P = P(1);
Cto = P / (8.314 * 10^-5) / T; 
E1 = 15000; 
E2 = 17500;
Ft = Fa + Fb + Fc + Fd; 
Ca = Cto * Fa / Ft;
Cb = Cto * Fb/Ft;
Cc = Cto * Fc/Ft;
Cd = Cto * Fd/Ft;
k1 = 0.075 * exp(E1 /1.987 * (1/300 - (1/T)));
k2 = 0.0015 * exp(E2 / 1.987 * (1 / 300 - (1 / T))); 
Fao = 50;  
rA = -(rA1 + rA2)
rB = -(rB1 + rB2 + rB3)
rC = rC1 + rC3
rD = rD1 + rD2 + rD3
rE = rE2 + rE3
% Differential equations
d(FA)/d(V) = rA
d(FB)/d(V) = rB
d(FC)/d(V) = rC
d(FD)/d(V) = rD
d(FE)/d(V) = rE
d(FF)/d(V) = rF
dYdV = [d(FC)/d(V)];
end

I didnt fully modify the script for your problem, but its something like that

[Torsten]

@Sean Powers:

And did the exaustive search lead to a different result than the usual optimization with "fminsearch" ?

[Bosnian Kingdom]

Thank you, I'll try.