Maximizing Batch Reactor Products in matlab

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]

Answer 1

[Torsten]

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

[naiva saeedia]

Tnx a lot u really saved my life and it has been pleasure to learn from u..I just have a question what if we change constraint to W_Ci-W_C10<0 instead of W_Ci-W_C10<=0.Since all W_Cis should be smaller than W_C10. Doesn't this change make code more accurate or it makes code to fail to solve the problem?

[Torsten]

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).

[naiva saeedia]

Tnx for ur answer. I understood ur point.

Answer 2

[Torsten]

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

[naiva saeedia]

Tnx a lot but actually I want to solve this with fmincon in form of optimization problem..I guess I wrote my question in a bad way..I want to change T,C_M,C_TEA and C_CAT at the same time to see how it effects x(1) to x(21) and then maximize x(15) compare to x(1),x(2),...x(14),x(16),...,x(21)..sorry for the miswritten question..can u please help me with this?

[Torsten]

Obviously, you have to choose T and C_TEA as small as possible and C_M and C_CAT as large as possible. No need to use "fmincon" for this.

[Torsten]

lb = [338.15;1;1;17];
ub = [368.15;100;200;400];
p0 = (lb+ub)/2;
p = fmincon(@fun,p0,[],[],[],[],lb,ub);

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

Topt = p(1)

Topt = 338.1500

C_M_opt = p(2)

C_M_opt = 100.0000

C_TEA_opt = p(3)

C_TEA_opt = 1.0000

C_CAT_opt = p(4)

C_CAT_opt = 400.0000

function value = fun(p)
  T=p(1);
  C_M=p(2);
  C_TEA=p(3);
  C_CAT=p(4);
  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]);
  value = -x(end,15);
end

[naiva saeedia]

So I tried what u said by changing T value to 338.15 and C_TEA to 1, C_M to 100 and C_CAT to 400..I got these results for x(3600,14), x(3600,15) and x(3600,16):

as u can see x(3600,15) is not greater than x(3600,14) so I guess I should use fmincon in order to maximize x(3600,15) and actually my professor insist that I should use MATLAB optimization library in order to solve this problem so If u could help me with this, It would be greatly appreciated and tnx a lot for ur comments!

[Torsten]

and what I meant by maximizing x(3600,15) is its value is going to be greater than x(3600,1),x(3600,2),...,x(3600,14),x(3600,15),...x(3600,21).

So the absolute value for x(3600,15) does not matter - it should only be greater than the other values ?

Or do you want to maximize x(3600,15) under the constraints that it should be greater than the other values ?

[naiva saeedia]

Yes the absolute value of x(3600,15) doesn't matter and it should only be greater than other values and we don't have constraints for values..all I want is x(3600,15)>x(3600,i) and i=1,2,..,14,16,..21...we have constraints for T,C_TEA,C_M and C_CAT...Here is a plot that I manage to find x(3600,15) max by changing C_TEA,C_M and C_CAT and T by hand....I performed some changes to x(3600,15) so the absolute max value is not the one that is here...Oh and one thing about C_CAT's constraint is it's C_CAT=x*C_TEA where x can change from 17 to 45..so C_TEA is not equal to 14 to 400(I typed that constraint wrongly sorry)...we have these constraint in general:

338.15=<T=<368.15

1=<C_M=<100

1=<C_TEA=<200

C_CAT=x*C_TEA where x=17,18,..,45

[Torsten]

Seems fmincon cannot find a combination of the 4 parameters that gives x(3600,15)>x(3600,i) and i=1,2,..,14,16,..21:

lb = [338.15;1;1;17];
ub = [368.15;100;200;400];
p0 = (lb+ub)/2;
options = optimoptions('fmincon','Display','iter','Algorithm','sqp');
p = fmincon(@obj,p0,[],[],[],[],lb,ub,@constr,options);

 Iter  Func-count            Fval   Feasibility   Step Length       Norm of   First-order  
                                                                       step    optimality
    0           5    0.000000e+00     9.310e+02     1.000e+00     0.000e+00     0.000e+00  
    1          10    0.000000e+00     3.155e+02     1.000e+00     9.145e+01     3.146e+02  
    2          15    0.000000e+00     4.830e+01     1.000e+00     1.181e+02     1.159e+02  
    3          20    0.000000e+00     3.069e+01     1.000e+00     2.472e+01     1.152e+02  
    4          25    0.000000e+00     3.093e+01     1.000e+00     2.054e+00     1.151e+02  
    5          30    0.000000e+00     3.092e+01     1.000e+00     7.485e-02     1.151e+02  
    6          31    0.000000e+00     3.092e+01     7.000e-01     7.397e-05     1.151e+02  

Converged to an infeasible point.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance but constraints are not
satisfied to within the value of the constraint tolerance.

Topt = p(1)

Topt = 368.1500

C_M_opt = p(2)

C_M_opt = 100

C_TEA_opt = p(3)

C_TEA_opt = 74.8944

C_CAT_opt = p(4)

C_CAT_opt = 17

function value = obj(p)
  %[t,x] = ode_solution(p);
  %value = -x(end,15);
  value = 0.0;
end

function [c,ceq] = constr(p)
  [t,x] = ode_solution(p);
  c = x(end,1:21)-x(end,15);
  ceq = [];
end

function [t,x] = ode_solution(p)
  T=p(1);
  C_M=p(2);
  C_TEA=p(3);
  C_CAT=p(4);
  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]);
end

[naiva saeedia]

So it doesn't have answer? right? tnx for ur help

[Torsten]

Would it be physically plausible that there exists such a combination ?

[naiva saeedia]

Actually I realized I have two mistakes..the objective funcion is not x(3600,15) but something else..since x(3600,15) cannot be maximize with these constraint..This is my code

clc
clear all
close all
T=  338.1500;
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= 100;
C_TEA=1;
C_CAT=400;
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;

%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;

actually the objective function is W_C10 and it should be greater than W_Ci where i=2,4,6,..,20 the constraint we have for this are listed below:

338.15=<T=<368.15

1=<C_M=<100

1=<C_TEA=<200

C_CAT=x*C_TEA where x=17,18,..,45

I tested this by hand and this objective function has a maximum with these constraints..I'm sure...the values that I found by hand are:

C_M=58

T=368.15

C_TEA=101

and

C_CAT=31*C_TEA

I'm so sorry for this problem..It's now compeletly correct

[Torsten]

W_C2 does not exist in your coding.

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','sqp');
  p = fmincon(@(p)obj(p,x(i),1),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);
  value_opt(i) = obj(p,x(i),2);
end
T_opt

T_opt = 1×29

364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524 364.2524

C_M_opt

C_M_opt = 1×29

46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898 46.8898

C_TEA_opt

C_TEA_opt = 1×29

102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056 102.3056

-value_opt

ans = 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

function value = obj(p,factor,iflag)
  [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;
  %if iflag == 1
  %  value = 0.0;
  %else
    value = -W_C10;
  %end
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]);
end

[naiva saeedia]

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]

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!