Maximizing Batch Reactor Products in matlab

4 次查看(过去 30 天)
Hi..I want to maximize W_C10 which is the mass concentraion of C10 products in a batch Reactor...And what I meant by maximization is W_C10 is going to be greater than all the other W_Ci,i=4,6,..,20..Here is my constraints for T,C_TEA and C_CAT and C_M:
338.15=<T=<368.15
1=<C_TEA<=200
1=<C_M<=100
C_CAT=x*C_TEA where x=17,18,...,45
Here is my ode system to solve W_Cis..Can someone help me to do this?
clc
clear
close all
T= 3.642755276270583e2;
T_r=273.15+0;
A1=4.525e5;
E1 =3985.;
A2 = 0.5101e9;
E2 = 52670.;
A3 = 0.4949e9;
E3 = 15140.;
A4 = 465000.;
E4 = 14330.;
A5 = 1.255;
E5 = 0;
R = 1.98;
k1 = A1*exp(E1*(1/T - 1/T_r)/R);
k2 = A2*exp(E2*(1/T - 1/T_r)/R);
k3 = A3*exp(E3*(1/T - 1/T_r)/R);
k4 = A4*exp(E4*(1/T - 1/T_r)/R);
k5 = A5*exp(E5*(1/T - 1/T_r)/R);
KA = 481.2;
KB = 277.1;
KC = 6906;
KD = 40830;
KE = 59.55;
C_M= 46.946030544906776;
C_TEA=1.022918116483290e2;
C_CAT=45*C_TEA;
C_M_K=C_M/(KB*C_TEA+1);
C_CAT_K=C_CAT/(KA*C_M+1);
C_TEA_k=C_TEA/(KC*C_CAT+1)^2;
dCdt=@(t,x) [k1.*C_CAT_K.*C_M_K-k2.*x(1).*C_M-k5.*x(1);k2.*C_M.*x(1)-(k3.*x(2).*C_M)./(KD.*C_TEA+1)-k5.*x(2);(k3.*x(2).*C_M)./(KD.*C_TEA+1)-(k3.*x(3).*C_M)./(KD.*C_TEA+1)-k5.*x(3);(k3.*x(3).*C_M)./(KD.*C_TEA+1)-(k3.*x(4).*C_M)./(KD.*C_TEA+1)-k5.*x(4);(k3.*x(4).*C_M)./(KD.*C_TEA+1)-(k3.*x(5).*C_M)./(KD.*C_TEA+1)-k5.*x(5);(k3.*x(5).*C_M)./(KD.*C_TEA+1)-(k3.*x(6).*C_M)./(KD.*C_TEA+1)-k5.*x(6);(k3.*x(6).*C_M)./(KD.*C_TEA+1)-(k3.*x(7).*C_M)./(KD.*C_TEA+1)-k5.*x(7);(k3.*x(7).*C_M)./(KD.*C_TEA+1)-(k3.*x(8).*C_M)./(KD.*C_TEA+1)-k5.*x(8);(k3.*x(8).*C_M)./(KD.*C_TEA+1)-(k3.*x(9).*C_M)./(KD.*C_TEA+1)-k5.*x(9);(k3.*x(9).*C_M)./(KD.*C_TEA+1)-(k3.*x(10).*C_M)./(KD.*C_TEA+1)-k5.*x(10);(k3.*x(10).*C_M)./(KD.*C_TEA+1)-(k3.*x(11).*C_M)./(KD.*C_TEA+1)-k5.*x(11);(k4.*x(3).*C_M)./(KE.*C_TEA+1)+k5.*x(3);(k4.*x(4).*C_M)./(KE.*C_TEA+1)+k5.*x(4);(k4.*x(5).*C_M)./(KE.*C_TEA+1)+k5.*x(5);(k4.*x(6).*C_M)./(KE.*C_TEA+1)+k5.*x(6);(k4.*x(7).*C_M)./(KE.*C_TEA+1)+k5.*x(7);(k4.*x(8).*C_M)./(KE.*C_TEA+1)+k5.*x(8);(k4.*x(9).*C_M)./(KE.*C_TEA+1)+k5.*x(9);(k4.*x(10).*C_M)./(KE.*C_TEA+1)+k5.*x(10);(k4.*x(11).*C_M)./(KE.*C_TEA+1)+k5.*x(11);k5.*(x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)+x(9)+x(10)+x(11))];
[t,x]=ode45(dCdt,[0,3600],[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0]);
%Moles Disturbutions in dead chains
t1=3600;
M_total=x(t1,12)+x(t1,13)+x(t1,14)+x(t1,15)+x(t1,16)+x(t1,17)+x(t1,18)+x(t1,19)+x(t1,20);
M_C4=x(t1,12)/M_total;
M_C6=x(t1,13)/M_total;
M_C8=x(t1,14)/M_total;
M_C10=x(t1,15)/M_total;
M_C12=x(t1,16)/M_total;
M_C14=x(t1,17)/M_total;
M_C16=x(t1,18)/M_total;
M_C18=x(t1,19)/M_total;
M_C20=x(t1,20)/M_total;
T1=table(M_C4,M_C6,M_C8,M_C10,M_C12,M_C14,M_C16,M_C18,M_C20)
%Mass Disturbutions in dead chains Basis=100mole
W_total=100*M_C4*72+100*M_C6*84+100*M_C8*112+100*M_C10*140+100*M_C12*168+100*M_C14*196+100*M_C16*224+100*M_C18*252+100*M_C20*280;
W_C4=100*M_C4*72/W_total;
W_C6=100*M_C6*84/W_total;
W_C8=100*M_C8*112/W_total;
W_C10=100*M_C10*140/W_total;
W_C12=100*M_C12*168/W_total;
W_C14=100*M_C14*196/W_total;
W_C16=100*M_C16*224/W_total;
W_C18=100*M_C18*252/W_total;
W_C20=100*M_C20*280/W_total;
T2=table(W_C4,W_C6,W_C8,W_C10,W_C12,W_C14,W_C16,W_C18,W_C20)
a=[W_C4,W_C6,W_C8,W_C10,W_C12,W_C14,W_C16,W_C18,W_C20]

采纳的回答

Torsten
Torsten 2023-8-2
编辑:Torsten 2023-8-2
You get almost identical values for the optimal T, C_M and C_TEA for all x-values in between 17 and 45 as you can see after the first loop in the code.
By strengthening the tolerances, now W_C4 is also smaller than W_C10 in all cases.
"obj" sets the function to be maximized - in this case I chose W_C10.
"constr" defines the constraints under which the function in "obj" is maximized. In this case, I chose W_Ci - W_C10 <= 0 for i = 4,6,...,20.
The code can be rearranged to look smarter (e.g. the threefold repetition of the calculation of the M_Ci ang W_Ci should be put in its own function) - I leave these technical changes to you.
x = 17:45;
for i = 1:numel(x)
lb = [338.15;1;1];
ub = [368.15;100;200];
p0 = (lb+ub)/2;
options = optimoptions('fmincon','Display','off','Algorithm','interior-point','ConstraintTolerance',1e-16);
[p,fval,exitflag,~] = fmincon(@(p)obj(p,x(i)),p0,[],[],[],[],lb,ub,@(p)constr(p,x(i)),options);
T_opt(i) = p(1);
C_M_opt(i) = p(2);
C_TEA_opt(i) = p(3);
C_CAT_opt(i) = C_TEA_opt(i)*x(i);
value_opt(i) = -fval;
exit(i) = exitflag;
end
exit
exit = 1×29
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
T_opt
T_opt = 1×29
364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755 364.2755
C_M_opt
C_M_opt = 1×29
46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460 46.9460
C_TEA_opt
C_TEA_opt = 1×29
102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918 102.2918
C_CAT_opt
C_CAT_opt = 1×29
1.0e+03 * 1.7390 1.8413 1.9435 2.0458 2.1481 2.2504 2.3527 2.4550 2.5573 2.6596 2.7619 2.8642 2.9665 3.0688 3.1710 3.2733 3.3756 3.4779 3.5802 3.6825 3.7848 3.8871 3.9894 4.0917 4.1940 4.2963 4.3985 4.5008 4.6031
value_opt
value_opt = 1×29
0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247 0.1247
%Postprocess results
for i = 1:numel(x)
[t,x] = ode_solution([T_opt(i),C_M_opt(i),C_TEA_opt(i)],x(i));
M_total=x(end,12)+x(end,13)+x(end,14)+x(end,15)+x(end,16)+x(end,17)+x(end,18)+x(end,19)+x(end,20);
M_C4=x(end,12)/M_total;
M_C6=x(end,13)/M_total;
M_C8=x(end,14)/M_total;
M_C10=x(end,15)/M_total;
M_C12=x(end,16)/M_total;
M_C14=x(end,17)/M_total;
M_C16=x(end,18)/M_total;
M_C18=x(end,19)/M_total;
M_C20=x(end,20)/M_total;
%Mass Disturbutions in dead chains Basis=100mole
W_total=100*M_C4*72+100*M_C6*84+100*M_C8*112+100*M_C10*140+100*M_C12*168+100*M_C14*196+100*M_C16*224+100*M_C18*252+100*M_C20*280;
W_C4(i)=100*M_C4*72/W_total;
W_C6(i)=100*M_C6*84/W_total;
W_C8(i)=100*M_C8*112/W_total;
W_C10(i)=100*M_C10*140/W_total;
W_C12(i)=100*M_C12*168/W_total;
W_C14(i)=100*M_C14*196/W_total;
W_C16(i)=100*M_C16*224/W_total;
W_C18(i)=100*M_C18*252/W_total;
W_C20(i)=100*M_C20*280/W_total;
end
any(W_C4>=W_C10)
ans = logical
0
any(W_C6>=W_C10)
ans = logical
0
any(W_C8>=W_C10)
ans = logical
0
any(W_C12>=W_C10)
ans = logical
0
any(W_C14>=W_C10)
ans = logical
0
any(W_C16>=W_C10)
ans = logical
0
any(W_C18>=W_C10)
ans = logical
0
any(W_C20>=W_C10)
ans = logical
0
function value = obj(p,factor)
[t,x] = ode_solution(p,factor);
M_total=x(end,12)+x(end,13)+x(end,14)+x(end,15)+x(end,16)+x(end,17)+x(end,18)+x(end,19)+x(end,20);
M_C4=x(end,12)/M_total;
M_C6=x(end,13)/M_total;
M_C8=x(end,14)/M_total;
M_C10=x(end,15)/M_total;
M_C12=x(end,16)/M_total;
M_C14=x(end,17)/M_total;
M_C16=x(end,18)/M_total;
M_C18=x(end,19)/M_total;
M_C20=x(end,20)/M_total;
%Mass Disturbutions in dead chains Basis=100mole
W_total=100*M_C4*72+100*M_C6*84+100*M_C8*112+100*M_C10*140+100*M_C12*168+100*M_C14*196+100*M_C16*224+100*M_C18*252+100*M_C20*280;
W_C4=100*M_C4*72/W_total;
W_C6=100*M_C6*84/W_total;
W_C8=100*M_C8*112/W_total;
W_C10=100*M_C10*140/W_total;
W_C12=100*M_C12*168/W_total;
W_C14=100*M_C14*196/W_total;
W_C16=100*M_C16*224/W_total;
W_C18=100*M_C18*252/W_total;
W_C20=100*M_C20*280/W_total;
value = -W_C10;
end
function [c,ceq] = constr(p,factor)
[t,x] = ode_solution(p,factor);
M_total=x(end,12)+x(end,13)+x(end,14)+x(end,15)+x(end,16)+x(end,17)+x(end,18)+x(end,19)+x(end,20);
M_C4=x(end,12)/M_total;
M_C6=x(end,13)/M_total;
M_C8=x(end,14)/M_total;
M_C10=x(end,15)/M_total;
M_C12=x(end,16)/M_total;
M_C14=x(end,17)/M_total;
M_C16=x(end,18)/M_total;
M_C18=x(end,19)/M_total;
M_C20=x(end,20)/M_total;
%Mass Disturbutions in dead chains Basis=100mole
W_total=100*M_C4*72+100*M_C6*84+100*M_C8*112+100*M_C10*140+100*M_C12*168+100*M_C14*196+100*M_C16*224+100*M_C18*252+100*M_C20*280;
W_C4=100*M_C4*72/W_total;
W_C6=100*M_C6*84/W_total;
W_C8=100*M_C8*112/W_total;
W_C10=100*M_C10*140/W_total;
W_C12=100*M_C12*168/W_total;
W_C14=100*M_C14*196/W_total;
W_C16=100*M_C16*224/W_total;
W_C18=100*M_C18*252/W_total;
W_C20=100*M_C20*280/W_total;
c(1) = W_C4 - W_C10;
c(2) = W_C6 - W_C10;
c(3) = W_C8 - W_C10;
c(4) = W_C12 - W_C10;
c(5) = W_C14 - W_C10;
c(6) = W_C16 - W_C10;
c(7) = W_C18 - W_C10;
c(8) = W_C20 - W_C10;
ceq = [];
end
function [t,x] = ode_solution(p,factor)
T=p(1);
C_M=p(2);
C_TEA=p(3);
C_CAT = factor*C_TEA;
T_r=273.15+0;
A1=4.525e5;
E1 =3985.;
A2 = 0.5101e9;
E2 = 52670.;
A3 = 0.4949e9;
E3 = 15140.;
A4 = 465000.;
E4 = 14330.;
A5 = 1.255;
E5 = 0;
R = 1.98;
k1 = A1*exp(E1*(1/T - 1/T_r)/R);
k2 = A2*exp(E2*(1/T - 1/T_r)/R);
k3 = A3*exp(E3*(1/T - 1/T_r)/R);
k4 = A4*exp(E4*(1/T - 1/T_r)/R);
k5 = A5*exp(E5*(1/T - 1/T_r)/R);
KA = 481.2;
KB = 277.1;
KC = 6906;
KD = 40830;
KE = 59.55;
C_M_K=C_M/(KB*C_TEA+1);
C_CAT_K=C_CAT/(KA*C_M+1);
C_TEA_k=C_TEA/(KC*C_CAT+1)^2;
dCdt=@(t,x) [k1.*C_CAT_K.*C_M_K-k2.*x(1).*C_M-k5.*x(1);k2.*C_M.*x(1)-(k3.*x(2).*C_M)./(KD.*C_TEA+1)-k5.*x(2);(k3.*x(2).*C_M)./(KD.*C_TEA+1)-(k3.*x(3).*C_M)./(KD.*C_TEA+1)-k5.*x(3);(k3.*x(3).*C_M)./(KD.*C_TEA+1)-(k3.*x(4).*C_M)./(KD.*C_TEA+1)-k5.*x(4);(k3.*x(4).*C_M)./(KD.*C_TEA+1)-(k3.*x(5).*C_M)./(KD.*C_TEA+1)-k5.*x(5);(k3.*x(5).*C_M)./(KD.*C_TEA+1)-(k3.*x(6).*C_M)./(KD.*C_TEA+1)-k5.*x(6);(k3.*x(6).*C_M)./(KD.*C_TEA+1)-(k3.*x(7).*C_M)./(KD.*C_TEA+1)-k5.*x(7);(k3.*x(7).*C_M)./(KD.*C_TEA+1)-(k3.*x(8).*C_M)./(KD.*C_TEA+1)-k5.*x(8);(k3.*x(8).*C_M)./(KD.*C_TEA+1)-(k3.*x(9).*C_M)./(KD.*C_TEA+1)-k5.*x(9);(k3.*x(9).*C_M)./(KD.*C_TEA+1)-(k3.*x(10).*C_M)./(KD.*C_TEA+1)-k5.*x(10);(k3.*x(10).*C_M)./(KD.*C_TEA+1)-(k3.*x(11).*C_M)./(KD.*C_TEA+1)-k5.*x(11);(k4.*x(3).*C_M)./(KE.*C_TEA+1)+k5.*x(3);(k4.*x(4).*C_M)./(KE.*C_TEA+1)+k5.*x(4);(k4.*x(5).*C_M)./(KE.*C_TEA+1)+k5.*x(5);(k4.*x(6).*C_M)./(KE.*C_TEA+1)+k5.*x(6);(k4.*x(7).*C_M)./(KE.*C_TEA+1)+k5.*x(7);(k4.*x(8).*C_M)./(KE.*C_TEA+1)+k5.*x(8);(k4.*x(9).*C_M)./(KE.*C_TEA+1)+k5.*x(9);(k4.*x(10).*C_M)./(KE.*C_TEA+1)+k5.*x(10);(k4.*x(11).*C_M)./(KE.*C_TEA+1)+k5.*x(11);k5.*(x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)+x(9)+x(10)+x(11))];
[t,x]=ode15s(dCdt,[0,3600],[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0],odeset('AbsTol',1e-10,'RelTol',1e-10));
end
  3 个评论
Torsten
Torsten 2023-8-2
编辑:Torsten 2023-8-2
It's not possible. Assume you want to minimize x^2 over x > 0. Then no optimal x exists. Only if you choose the constraint as x>=0, you get x=0 as solution.
This generalizes to the domains over which optimization is performed in general. It has to be always closed (<=), not open (<0).
But I changed the checks to any(W_Ci>=W_C10) in the above code, and as you can see: all constraints are satisfied in the strict sense (W_Ci < W_C10).

请先登录,再进行评论。

更多回答(1 个)

Torsten
Torsten 2023-8-1
编辑:Torsten 2023-8-1
max_value = -Inf;
T_array = linspace(338.15,368.15,10);
C_M_array = linspace(1,100,10);
C_TEA_array = linspace(1,200,10);
C_CAT_array = linspace(17,400,10);
for i = 1:numel(T_array)
T = T_array(i);
for j = 1:numel(C_M_array)
C_M = C_M_array(j);
for k = 1:numel(C_TEA_array)
C_TEA = C_TEA_array(k);
for l = 1:numel(C_CAT_array)
C_CAT = C_CAT_array(l);
T_r=273.15+0;
A1=4.525e5;
E1 =3985.;
A2 = 0.5101e9;
E2 = 52670.;
A3 = 0.4949e9;
E3 = 15140.;
A4 = 465000.;
E4 = 14330.;
A5 = 1.255;
E5 = 0;
R = 1.98;
k1 = A1*exp(E1*(1/T - 1/T_r)/R);
k2 = A2*exp(E2*(1/T - 1/T_r)/R);
k3 = A3*exp(E3*(1/T - 1/T_r)/R);
k4 = A4*exp(E4*(1/T - 1/T_r)/R);
k5 = A5*exp(E5*(1/T - 1/T_r)/R);
KA = 481.2;
KB = 277.1;
KC = 6906;
KD = 40830;
KE = 59.55;
C_M_K=C_M/(KB*C_TEA+1);
C_CAT_K=C_CAT/(KA*C_M+1);
C_TEA_k=C_TEA/(KC*C_CAT+1)^2;
dCdt=@(t,x) [k1.*C_CAT_K.*C_M_K-k2.*x(1).*C_M-k5.*x(1);k2.*C_M.*x(1)-(k3.*x(2).*C_M)./(KD.*C_TEA+1)-k5.*x(2);(k3.*x(2).*C_M)./(KD.*C_TEA+1)-(k3.*x(3).*C_M)./(KD.*C_TEA+1)-k5.*x(3);(k3.*x(3).*C_M)./(KD.*C_TEA+1)-(k3.*x(4).*C_M)./(KD.*C_TEA+1)-k5.*x(4);(k3.*x(4).*C_M)./(KD.*C_TEA+1)-(k3.*x(5).*C_M)./(KD.*C_TEA+1)-k5.*x(5);(k3.*x(5).*C_M)./(KD.*C_TEA+1)-(k3.*x(6).*C_M)./(KD.*C_TEA+1)-k5.*x(6);(k3.*x(6).*C_M)./(KD.*C_TEA+1)-(k3.*x(7).*C_M)./(KD.*C_TEA+1)-k5.*x(7);(k3.*x(7).*C_M)./(KD.*C_TEA+1)-(k3.*x(8).*C_M)./(KD.*C_TEA+1)-k5.*x(8);(k3.*x(8).*C_M)./(KD.*C_TEA+1)-(k3.*x(9).*C_M)./(KD.*C_TEA+1)-k5.*x(9);(k3.*x(9).*C_M)./(KD.*C_TEA+1)-(k3.*x(10).*C_M)./(KD.*C_TEA+1)-k5.*x(10);(k3.*x(10).*C_M)./(KD.*C_TEA+1)-(k3.*x(11).*C_M)./(KD.*C_TEA+1)-k5.*x(11);(k4.*x(3).*C_M)./(KE.*C_TEA+1)+k5.*x(3);(k4.*x(4).*C_M)./(KE.*C_TEA+1)+k5.*x(4);(k4.*x(5).*C_M)./(KE.*C_TEA+1)+k5.*x(5);(k4.*x(6).*C_M)./(KE.*C_TEA+1)+k5.*x(6);(k4.*x(7).*C_M)./(KE.*C_TEA+1)+k5.*x(7);(k4.*x(8).*C_M)./(KE.*C_TEA+1)+k5.*x(8);(k4.*x(9).*C_M)./(KE.*C_TEA+1)+k5.*x(9);(k4.*x(10).*C_M)./(KE.*C_TEA+1)+k5.*x(10);(k4.*x(11).*C_M)./(KE.*C_TEA+1)+k5.*x(11);k5.*(x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)+x(9)+x(10)+x(11))];
[t,x]=ode15s(dCdt,[0,3600],[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0]);
if x(end,15) > max_value
I = i;
J = j;
K = k;
L = l;
max_value = x(end,15);
end
end
end
end
end
T_array(I)
ans = 338.1500
C_M_array(J)
ans = 100
C_TEA_array(K)
ans = 1
C_CAT_array(L)
ans = 400
max_value
max_value = 9.9395e+05
  13 个评论
naiva saeedia
naiva saeedia 2023-8-2
Wow..tnx a lot and Yeah it didn't have WC2 sorry for the mistake...Can we find what's the value of factor or C_CAT is here?
naiva saeedia
naiva saeedia 2023-8-2
编辑:naiva saeedia 2023-8-2
Never mind...the value of factor is 29 right?...there is a little problem and I don't know why this is happening...I tried the value that your code is giving me in my code to find out WCi's values but as u can see WC4 is greater than WC10
Why this is happening? Every other WCis are lesser than the WC10 and WC4 is the only value that is greater than W_C10...I also didn't understand what Constr and Obj functions are doing so If u could give me some explanation that would be great(especially when u defined C(1),..,C(8) in Constr and value part in Obj function) and tnq for ur time and effort what u did is pretty amazing!

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Solver Outputs and Iterative Display 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by