How to apply Genetic Algorithm solver to do parameter estimation for ODE based model?

2 次查看（过去 30 天）
Bendo Voon
Bendo Voon 2014-7-12
评论： Jan 2019-3-14
I have a ODE based model and a set of experimental data for my project. There are seven parameters (Constant) inside the model need to do parameter estimation through Genetic Algorithm solver. My fitness function/Objective function will be sum of quadratic difference between simulated and experimental result.
pm(..) is the parameter I need to do estimation to best fit my curve. The function cgtase_cont is the mathematical model used. I tried to use GA solver to run it but it always appeared "Optimization running. Error running optimization. Too many output arguments.". Function square_error_ga is my objective function.
Using GA to solve this problem is the requirement of the project. Thanks in advance!
----------------------------------------------------------------------
Function cgtase_cont code:
function dx = cgtase_cont(t,x,ux,Ks,K1,H,KI,Kox,Ag,Ea,R,Ad,Ed,dH,dS,dCp,Temp0,rol,N,Di,V,Vs,...
K2,Kd,Ysx,mx,Clx,qO2,YH,Ah,Cp,U,uxp,I,Ki,K,D,Tempj)
dx = zeros(6,1);
%state variables
u = (ux*x(2))/((Ks+x(2))*(1+K1/H+H/K2)*(x(2)+KI+(x(2)^2)/Ks)*(Kox+1.4e-3))*(Ag*exp(-Ea/(R*x(4)))-Ad*exp(-Ed/(R*x(4))));
dGTS = dH-x(4)*dS+dCp*(x(4)-Temp0-x(4)*log(x(4)/Temp0));
Pg = 0.4*6*rol*N^3*Di^5;
Kla = 0.0195*((Pg/V)^0.55)*((Vs)^0.64)*(1+2.12*x(1)+0.2*x(1)^2)^-0.25;
Kn = 111.3*exp(-((x(4)-313.3)/7.42)^2);
%models
dx(1) = (u-Kd)*x(1);
dx(2) = (u/Ysx+mx)*(-x(1));
dx(3) = (Kla*(Clx-x(3))-qO2*x(1));
dx(4) = ((u*x(1))/(YH*rol*Cp)-(U*Ah)/(V*rol*Cp)*(x(4)-Tempj));
dx(5) = ((uxp*I)/(I+Ki)*u*x(1)-K*x(5));
dx(6) = (D*exp(-dGTS/(R*x(4))))*x(5)-Kn*x(6);
------------------------------------------------------------
Function square_error_ga code:
function error = squared_error_ga(data_time, data_points, fn_time, fn_points)
% Intepolate the fit to match the time points of the data
fn_values = interp1(fn_time, fn_points, data_time);
% Square the error values and sum them up
error = sum((fn_values - data_points).^2);
--------------------------------------------------------
Function ga_error code:
function error = ga_error(pm)
%Estimated Parameter set
Ks = pm(1); %half saturation constant (g/l) pm(1)
Ki = pm(2); %saturation constant for inducer (g/l) pm(2)
Ag=10^pm(3); %Ag = 10^7.9; %Arrhenius constant for growth pm(3)
Ad=10^pm(4); %Ad = 10^10; %Arrhenius constant for death pm(4)
KI= pm(5); %inhibition constant due to excessive subtrate (g/l)pm(5)
mx = pm(6); %cell maintenance coefficient on substrate pm(6)
Ysx = pm(7); %yield coefficient (g/g)pm(7)
D = pm(8); %preexponential factor (per h) pm(8)
Kox = pm(9); %oxygen limitation constant (g/l) pm(9)
%Actual Constant
%Ks = 7.5; %half saturation constant (g/l)
%Ki = 0.001; %saturation constant for inducer (g/l)
%Ag = 10^7.9; %Arrhenius constant for growth
%Ad = 10^10; %Arrhenius constant for death
%KI= 0.001; %inhibition constant due to excessive subtrate (g/l)
%mx = 0.003; %cell maintenance coefficient
%Ysx = 0.0105; %yield coefficient (g/g)
%D = 1.65*10^14; %preexponential factor (per h)
%Kox = 7.104*10^-7; %oxygen limitation constant (g/l
ux = 0.1027; %maximum specific growth rate (per h)
K1 = 7.05*10^-10; %constant (Molar)
K2 = 6.86*10^-8; %constant (Molar)
H = 10^(-7.1); %ion hydrogen concentration (M)
Ea = 10000; %activation energy for growth (cal/mol)
Ed = 80000; %activation energy for death (cal/mol)
R = 1.987; %gas constant (cal/mol K)
uxp = 0.5872; %maximum specific production rate (per h)
I = 0.162; %inducer concentration (g/l)
K = 0.0625; %protein degradation rate
dH = 9.4*10^3; %change in enthalpy (cal/mol)
dS = -0.01*10^3; %change in entropy (cal/mol)
dCp = -0.32*10^3; %change in the heat capacity of the protein between fold and unfold state
Temp0 = 295; %reference temperature (K)
rol = 998.23; %density medium (g/l)
Di = 46*10^-3; %impeller diameter (m), Rushton
V = 1.5; %volume culture (l)
Vs = 1.1372; %oxygen superficial velocity (m/s)
Clx = 7*10^-3; %saturated DO concentration (g/l)
qO2 = 384*10^-3; %specific uptake rate of oxygen (g/g h)
U = 3.6482*10^4; %overall heat transfer coefficient (cal/h m2 K)
Ah = 0.606; %contacted surface area (m2)
YH = 0.42*10^-3; %1/YH is metabolic heat evolved per g biomass (g/cal)
Cp = 1; %heat capacity of medium (cal/ g K)
N = 10/3; %agitation speed (rev/s)
Kd = 0.02; %specific death rate (per h)
Tempj=293;
%-------------------------------------------
ode45(@(t,x)cgtase_cont(t,x,ux,Ks,K1,H,KI,Kox,Ag,Ea,R,Ad,Ed,dH,dS,dCp,Temp0,rol,N,Di,V,Vs,...
K2,Kd,Ysx,mx,Clx,qO2,YH,Ah,Cp,U,uxp,I,Ki,K,D,Tempj), ...
[0 60], [0.12 50 1.4e-3 310 0 0]);
error = squared_error_ga([0 4 8 12 16 20 24 28 32 36 40 44 48] , ...
[0 0 0.7586 1.2944 4.9136 6.9158 10.3277...
10.0373 10.8607 12.9424 12.6819 11.7554 10.6795], ...
t, ...
y(:,6));
------------------------------------------------------
The graph of my project is something like:
%
  2 个评论
Jan
Jan 2019-3-14
@Joker: Really exactly the same question? In the original question, the output of ode45 was not caught, so I cannot guess, why the code runs at all. As far as I can see, the GA-part has not been posted and the error message was shown partially only. So please be so kind and open a new thread providing the required details to help you.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Genetic Algorithm 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by