Can genetic algorithm be used to find two independent optimum operating conditions for predefined input and output ?

I've built a moving bed biomass torrefaction reactor which is controlled by the speed of the conveying mechanism and the air mass flow rate fed for combustion with the released volatiles to give a certain temperature. I have included all the related equations of reactor behavior into a straightforward Matlab code where I define the input characteristics of biomass and the operating conditions and then the code gives the output characteristics for such conditions. Now I want to get the optimum operating conditions (time in terms of conveyor speed and temperature in terms of air mass flow rate) to reach a predefined output (biomass characteristics). I'm not familiar with optimization techniques but after some research I got to know that genetic algorithm is the most used algorithm for optimization. So, I'm asking if it can be of help for my concern and if so, is there a guide to do such optimization?

[Star Strider]

Optimisation functions, including the genetic algorithm (ga function) optimise a set of parameters to produce the desired result. It is certainly possible to create a fitness function that will optimise your model, however it would be necessary to see the model in order to determine the best approach.

[Mohamed Talaat]

just forget the details of the above question and let's focus on the code I have built itself. this is the code

clear;
t0=0;   %%%%% Starting time when Temp. is 200 deg C %%%%%%%%
T_initial=200;     %%%%%%% Initial Temp. %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% disp ('Please enter the following data ...');
% prompt = 'Ea1 = ';      Ea1 = input(prompt);      %%%% Activation energy is denoted Ea %%%%%%%%
% prompt = 'K01 = ';      K01 = input(prompt);      %%%% The pre-exponential factor is denoted K0 %%%%%%%
% prompt = 'EaV1 = ';     EaV1 = input(prompt);
% prompt = 'K0V1 = ';     K0V1 = input(prompt);
% prompt = 'Ea2 = ';      Ea2 = input(prompt);
% prompt = 'K02 = ';      K02 = input(prompt);
% prompt = 'EaV2 = ';     EaV2 = input(prompt);
% prompt = 'K0V2 = ';     K0V2 = input(prompt);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Constants of Willow reactions %%%%%%%%
Ea1 = 75976;
K01 = 24800;
EaV1 = 114214;
K0V1 = 32300000;
Ea2 = 151711;
K02 = 11000000000;
EaV2 = 151711;
K0V2 = 15900000000;
%%%%%%%                      %%%%%%
%%%%%%% V1 ultimate analysis %%%%%%
Yv1_c=0.181155576;       
Yv1_H=0.072600424;
Yv1_o=0.746243999;
Yv1_V=0;
Yv1_ash=0;
%%%%%%%                      %%%%%%
%%%%%%% V2 ultimate analysis %%%%%%
Yv2_c=0.343756522;
Yv2_H=0.088292754;
Yv2_o=0.568950725;
Yv2_N=0;
Yv2_ash=0;
%%%%%%%                      %%%%%%
%%%%%%%  A ultimate analysis %%%%%%
YA_c=0.472;
YA_H=0.061;
YA_o=0.451;
YA_N=0.003;
YA_S=0;
YA_ash=0.013;
%%%%%%%                      %%%%%%
YB_ash=0;
YB_N=0;
YB_S=0;

YC_ash=0;
YC_N=0;
YC_S=0;

%%%%%%%    A (C_a H_b O_c)   %%%%%%
A_C_a=YA_c/0.012011;
A_H_b=YA_H/0.001079;
A_O_c=YA_o/0.016;
%%%%%%%                      %%%%%%
%%%%%%%                      %%%%%%
%%%%%%% Complete combustion of A gives the following equation:
%%%%%%% Ca Hb Oc + x O2 = y H2O + z CO2 %%%%%%
%%%%%%% x = (2a + b/2 - c)/2
%%%%%%% y = b/2
%%%%%%% z = a
%%%%%%%                      %%%%%%
A_x=(2*A_C_a+A_H_b/2-A_O_c)/2;
A_y=A_H_b/2;
A_z=A_C_a;

HHV_A_Boie = 1000*[351.69*(YA_c) + 1162.46*(YA_H) - 110.95*(YA_o) + 104.67*(YA_S) - 62.8*(YA_N)] ;
HHV_A_Friedl = 10^6*[35.45*(YA_c)^2 - 23.15*(YA_c) - 223.35*(YA_H) + 512*(YA_c)*(YA_H) + 13.1*(YA_N) + 20.5875] ;
HHV_A_IGT = 1000*[341.7*(YA_c) + 1322.1*(YA_H) - 119.8*(YA_o + YA_N) - 15.3*(YA_ash)] ;

prompt = 'Heating rate (degrees/min) = ';     deg_min = input(prompt);
prompt = 'Final temperature = ';     T_final = input(prompt);
disp ('Heating up period is '); disp ((T_final-T_initial)/deg_min);
prompt = 'Isothermal heating period (min) = ';     t_iso = input(prompt);
prompt = 'Solution time step (sec) = ';     dt = input(prompt);
% t_total=(T_final-T_initial)/deg_min+t_iso;
dT=deg_min/60*dt;      %%%%% Temp. difference in one time step %%%%%%%
N=((T_final-T_initial)/deg_min+t_iso)*60/dt;    %%%%%%%%%%% Number of time steps %%%%%%%%%%%
A= ones(N,1);
B= zeros(N,1);
C= zeros(N,1);
V1= zeros(N,1);
V2= zeros(N,1);
Ysolid_c=YA_c*ones(N,1);
Ysolid_H=YA_H*ones(N,1);
Ysolid_o=YA_o*ones(N,1);
Ysolid_ash=YA_ash*ones(N,1);
RHHV_Boie=ones(N,1);
RHHV_Friedl=ones(N,1);
RHHV_IGT=ones(N,1);
Va=zeros(N,1);
Vb=zeros(N,1);
Vc=zeros(N,1);
Vd=zeros(N,1);
Ve=zeros(N,1);
Vf=zeros(N,1);
Vg=zeros(N,1);
Vh=zeros(N,1);
Vi=zeros(N,1);
count =-1;
for i = 2 : N+1
    count =  count+1;
    T(i-1,1) = T_initial + count*dT;
    if (T(i-1,1)>T_final)
        T(i-1,1)=T_final;
    end
    K1(i-1,1)= K01*exp(-Ea1/8.314/(T(i-1)+273));
    KV1(i-1,1)= K0V1*exp(-EaV1/8.314/(T(i-1)+273));
    K2(i-1,1)= K02*exp(-Ea2/8.314/(T(i-1)+273));
    KV2(i-1,1)= K0V2*exp(-EaV2/8.314/(T(i-1)+273));
    ra(i-1,1)=-(K1(i-1)+KV1(i-1))*A(i-1);
    rb(i-1,1)=K1(i-1)*A(i-1)-(K2(i-1)+KV2(i-1))*B(i-1);
    rc(i-1,1)=K2(i-1)*B(i-1);
    rv1(i-1,1)=KV1(i-1)*A(i-1);
    rv2(i-1,1)=KV2(i-1)*B(i-1);
    beta(i-1,1)=K1(i-1,1)/(K1(i-1,1)+KV1(i-1,1));
    nu(i-1,1)=KV1(i-1,1)/(K1(i-1,1)+KV1(i-1,1));
    gamma(i-1,1)=K2(i-1,1)/(K2(i-1,1)+KV2(i-1,1));
    zeta(i-1,1)=KV2(i-1,1)/(K2(i-1,1)+KV2(i-1,1));
    YB_c(i-1,1)=(YA_c-nu(i-1,1)*Yv1_c)/beta(i-1,1);
    YB_H(i-1,1)=(YA_H-nu(i-1,1)*Yv1_H)/beta(i-1,1);
    YB_o(i-1,1)=(YA_o-nu(i-1,1)*Yv1_o)/beta(i-1,1);
    YC_c(i-1,1)=(YB_c(i-1,1)-zeta(i-1,1)*Yv2_c)/gamma(i-1,1);
    YC_H(i-1,1)=(YB_H(i-1,1)-zeta(i-1,1)*Yv2_H)/gamma(i-1,1);
    YC_o(i-1,1)=(YB_o(i-1,1)-zeta(i-1,1)*Yv2_o)/gamma(i-1,1);
    Ysolid_c(i,1)=Ysolid_c(i-1,1)-dt*(Yv1_c*rv1(i-1,1)+Yv2_c*rv2(i-1,1));
    Ysolid_H(i,1)=Ysolid_H(i-1,1)-dt*(Yv1_H*rv1(i-1,1)+Yv2_H*rv2(i-1,1));
    Ysolid_o(i,1)=Ysolid_o(i-1,1)-dt*(Yv1_o*rv1(i-1,1)+Yv2_o*rv2(i-1,1));
    Ysolid_ash(i,1)=Ysolid_ash(i-1,1)-dt*(Yv1_ash*rv1(i-1,1)+Yv2_ash*rv2(i-1,1));
    
    HHV_B_Boie(i-1,1) = 1000*[351.69*(YB_c(i-1,1)) + 1162.46*(YB_H(i-1,1)) - 110.95*(YB_o(i-1,1)) + 104.67*(YB_S) - 62.8*(YB_N)] ;
    HHV_B_Friedl(i-1,1) = 10^6*[35.45*(YB_c(i-1,1))^2 - 23.15*(YB_c(i-1,1)) - 223.35*(YB_H(i-1,1)) + 512*(YB_c(i-1,1))*(YB_H(i-1,1)) + 13.1*(YB_N) + 20.5875] ;
    HHV_B_IGT(i-1,1) = 1000*[341.7*(YB_c(i-1,1)) + 1322.1*(YB_H(i-1,1)) - 119.8*(YB_o(i-1,1) + YB_N) - 15.3*(YB_ash)] ;
    HHV_C_Boie(i-1,1) = 1000*[351.69*(YC_c(i-1,1)) + 1162.46*(YC_H(i-1,1)) - 110.95*(YC_o(i-1,1)) + 104.67*(YC_S) - 62.8*(YC_N)] ;
    HHV_C_Friedl(i-1,1) = 10^6*[35.45*(YC_c(i-1,1))^2 - 23.15*(YC_c(i-1,1)) - 223.35*(YC_H(i-1,1)) + 512*(YC_c(i-1,1))*(YC_H(i-1,1)) + 13.1*(YC_N) + 20.5875] ;
    HHV_C_IGT(i-1,1) = 1000*[341.7*(YC_c(i-1,1)) + 1322.1*(YC_H(i-1,1)) - 119.8*(YC_o(i-1,1) + YC_N) - 15.3*(YC_ash)] ;
        
    Va(i,1)=Va(i-1,1)+dt*(0.148*rv1(i-1,1)+0.161*rv2(i-1,1));
    Vb(i,1)=Vb(i-1,1)+dt*(0.481*rv1(i-1,1)+0.076*rv2(i-1,1));
    Vc(i,1)=Vc(i-1,1)+dt*(0.053*rv1(i-1,1)+0.051*rv2(i-1,1));
    Vd(i,1)=Vd(i-1,1)+dt*(0.042*rv1(i-1,1)+0.301*rv2(i-1,1));
    Ve(i,1)=Ve(i-1,1)+dt*(0.013*rv1(i-1,1)+0.313*rv2(i-1,1));
    Vf(i,1)=Vf(i-1,1)+dt*(0.011*rv1(i-1,1)+0.000*rv2(i-1,1));
    Vg(i,1)=Vg(i-1,1)+dt*(0.006*rv1(i-1,1)+0.097*rv2(i-1,1));
    Vh(i,1)=Vh(i-1,1)+dt*(0.204*rv1(i-1,1)+0.000*rv2(i-1,1));
    Vi(i,1)=Vi(i-1,1)+dt*(0.042*rv1(i-1,1)+0.001*rv2(i-1,1));
 
    dVa_dt(i-1,1)=60*(0.148*rv1(i-1,1)+0.161*rv2(i-1,1));
    dVb_dt(i-1,1)=60*(0.481*rv1(i-1,1)+0.076*rv2(i-1,1));
    dVc_dt(i-1,1)=60*(0.053*rv1(i-1,1)+0.051*rv2(i-1,1));
    dVd_dt(i-1,1)=60*(0.042*rv1(i-1,1)+0.301*rv2(i-1,1));
    dVe_dt(i-1,1)=60*(0.013*rv1(i-1,1)+0.313*rv2(i-1,1));
    dVf_dt(i-1,1)=60*(0.011*rv1(i-1,1)+0.000*rv2(i-1,1));
    dVg_dt(i-1,1)=60*(0.006*rv1(i-1,1)+0.097*rv2(i-1,1));
    dVh_dt(i-1,1)=60*(0.204*rv1(i-1,1)+0.000*rv2(i-1,1));
    dVi_dt(i-1,1)=60*(0.042*rv1(i-1,1)+0.001*rv2(i-1,1));
    
    A(i,1)=A(i-1,1)-dt*(A(i-1,1)*(K1(i-1,1)+KV1(i-1,1)));
    B(i,1)=B(i-1,1)+dt*(A(i-1,1)*K1(i-1,1)-B(i-1,1)*(K2(i-1,1)+KV2(i-1,1)));
    C(i,1)=C(i-1,1)+dt*B(i-1,1)*K2(i-1,1);
    V1(i,1)=V1(i-1,1)+dt*A(i-1,1)*KV1(i-1,1);
    V2(i,1)=V2(i-1,1)+dt*B(i-1,1)*KV2(i-1,1);
    t_total(i,1)=i*dt/60;
    V_yield_cond=Va+Vb+Vc+Vd+Ve+Vf+Vg;
    V_yield_non_cond=Vh+Vi;
    Solid_loss=1-(A+B+C);
    
    RHHV_Boie(i,1)=(A(i,1)*HHV_A_Boie + B(i,1)*HHV_B_Boie(i-1,1) + C(i,1)*HHV_C_Boie(i-1,1))/((A(i,1)+B(i,1)+C(i,1))*HHV_A_Boie) ;
    RHHV_Friedl(i,1)=(A(i,1)*HHV_A_Friedl + B(i,1)*HHV_B_Friedl(i-1,1) + C(i,1)*HHV_C_Friedl(i-1,1))/((A(i,1)+B(i,1)+C(i,1))*HHV_A_Friedl) ;
    RHHV_IGT(i,1)=(A(i,1)*HHV_A_IGT + B(i,1)*HHV_B_IGT(i-1,1) + C(i,1)*HHV_C_IGT(i-1,1))/((A(i,1)+B(i,1)+C(i,1))*HHV_A_IGT) ;
    
end
A=A([1:N],:);
B=B([1:N],:);
C=C([1:N],:);
V1=V1([1:N],:);
V2=V2([1:N],:);
ra=ra([1:N],:);
rb=rb([1:N],:);
rc=rc([1:N],:);
rv1=rv1([1:N],:);
rv2=rv2([1:N],:);
Ysolid_c=Ysolid_c([1:N],:);
Ysolid_H=Ysolid_H([1:N],:);
Ysolid_o=Ysolid_o([1:N],:);
Ysolid_ash=Ysolid_ash([1:N],:);
Va=Va([1:N],:);
Vb=Vb([1:N],:);
Vc=Vc([1:N],:);
Vd=Vd([1:N],:);
Ve=Ve([1:N],:);
Vf=Vf([1:N],:);
Vg=Vg([1:N],:);
Vh=Vh([1:N],:);
Vi=Vi([1:N],:);
RHHV_Boie=RHHV_Boie([1:N],:);
RHHV_Friedl=RHHV_Friedl([1:N],:);
RHHV_IGT=RHHV_IGT([1:N],:);

V_yield_cond=V_yield_cond([1:N],:);
V_yield_non_cond=V_yield_non_cond([1:N],:);
Solid_loss=Solid_loss([1:N],:);
Ysolid_c_perc=Ysolid_c/YA_c;
Ysolid_H_perc=Ysolid_H/YA_H;
Ysolid_o_perc=Ysolid_o/YA_o;
YS=1-Solid_loss;

t_total=t_total([1:N],:);
comb = [];
for k1 = 1:1
    comb = [comb Ysolid_c(:,k1) Ysolid_H(:,k1) Ysolid_o(:,k1) Ysolid_ash(:,k1)];
end
comb2 = [];
for k2 = 1:1
    comb2 = [comb2 Va(:,k2) Vb(:,k2) Vc(:,k2) Vd(:,k2) Ve(:,k2) Vf(:,k2) Vg(:,k2) Vh(:,k2) Vi(:,k2) ];
end

% figure(1);
% area (t_total,comb);
% figure(2);
% area (t_total,comb2);
% figure(3);
% plot(t_total,dVa_dt,t_total,dVb_dt,t_total,dVc_dt,t_total,dVd_dt,t_total,dVe_dt,t_total,dVf_dt,t_total,dVg_dt,t_total,dVh_dt,t_total,dVi_dt);
% figure(4);
% plot(Solid_loss,V_yield_cond,Solid_loss,V_yield_non_cond);
% figure(5);
% plot(T,beta,T,nu,T,gamma,T,zeta);
% figure(6);
% plot(T,YB_c,T,YB_H,T,YB_o);
% figure(7);
% plot(T,YC_c,T,YC_H,T,YC_o);
% figure(8);
% plot(Solid_loss,Ysolid_c_perc,Solid_loss,Ysolid_H_perc,Solid_loss,Ysolid_o_perc);
% figure(9);
% plot(Solid_loss,RHHV_Boie,'--',Solid_loss,RHHV_Friedl);

I enter the heating up rate in terms of deg/min (~10) and the final temperature (~200-300) and the isothermal time (~10-30) and the time step size and the code runs and gives me the results in this case (I mean these times and temperatures). Now I want to get the time (total) and final temperature in case I need a specific results in terms of (Ysolid_c & Ysolid_H & Ysolid_o)

Answer 1

[Alan Weiss]

It sounds to me as if you are trying to optimize an ODE system, possibly fitting the parameters to an existing function. If I understand that correctly, then I suggest that you look at these examples:

Fit an Ordinary Differential Equation (ODE)

Fit ODE, Problem-Based

Optimize an ODE in Parallel

Alan Weiss

MATLAB mathematical toolbox documentation

[Mohamed Talaat]

Dear Alan, you can find the code in comments, so lock at it and kindly give your suggestions