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

1 ビュー (過去 30 日間)
Mohamed Talaat 2019 年 10 月 29 日
Commented: Mohamed Talaat 2019 年 11 月 1 日
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 件のコメント

Star Strider 2019 年 10 月 29 日
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 2019 年 11 月 1 日
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)

サインイン to comment.

### 件の回答 (1)

Alan Weiss 2019 年 10 月 29 日
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 件のコメント

Mohamed Talaat 2019 年 11 月 1 日
Dear Alan, you can find the code in comments, so lock at it and kindly give your suggestions

サインイン to comment.

サインイン してこの質問に回答します。