MATLAB Answers

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

2 views (last 30 days)
Mohamed Talaat
Mohamed Talaat on 29 Oct 2019
Commented: Mohamed Talaat on 1 Nov 2019
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?

  2 Comments

Star Strider
Star Strider on 29 Oct 2019
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
Mohamed Talaat on 1 Nov 2019
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)

Sign in to comment.

Answers (1)

Alan Weiss
Alan Weiss on 29 Oct 2019
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:
Alan Weiss
MATLAB mathematical toolbox documentation

  1 Comment

Sign in to comment.

Sign in to answer this question.


Translated by