Incorrect outcome by applying Genetic Algorithm solver to do parameter estimation for ODE based model

1 次查看（过去 30 天）
Bendo Voon
Bendo Voon 2014-7-12
编辑： Bendo Voon 2014-7-12
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_model is the mathematical model used. Function square_error_ga is my objective function.I tried to use GA solver to run it and it can successfully run it but the parameter value is totally different with the original or guess value. Besides, the fitness value chart is always a straight line showed as below.
Can anyone point out my mistakes on my coding? Thanks in advance.
function error = ga_error(pm)
load enzyme_bendo.mat
global um K1 K2 H Kox KI Ag Ad Ea Ed R ug ug1 ug2 ug3 X uxp I K D dH...
dS dCp Temp0 dGTS Pn Pc Ysx mx S rol Di V Vs Clx qO2 Pg Kla Cl Fj Vj...
Tempj0 U Ah rolj Cpj Tempj YH Cp Temp Kd Kn t N
%Actual Constant
%Ks = 7.5; %half saturation constant (g/l) pm(1)
%Ki = 0.001; %saturation constant for inducer (g/l) pm(2)
%Ag = 10^7.9; %Arrhenius constant for growth pm(3)
%Ad = 10^10; %Arrhenius constant for death pm(4)
%KI= 0.001; %inhibition constant due to excessive subtrate (g/l) pm(5)
%mx = 0.003; %cell maintenance coefficient pm(6)
%Ysx = 0.0105; %yield coefficient (g/g) pm(7)
%D = 1.65*10^14; %preexponential factor (per h) pm(8)
%Kox = 7.104*10^-7; %oxygen limitation constant (g/l) pm(9)
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;
% Call to odefun
%options = odeset('RelTol', 1e-4,'NonNegative', [1 2 3]);
[t,y] = ode45(@(t,x)cgtase_model(t,x,ux,K1,H,Ea,R,Ed,dH,dS,dCp,Temp0,rol,N,Di,V,Vs,...
K2,Kd,Clx,qO2,YH,Ah,Cp,U,uxp,I,K,Tempj,pm), ...
[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));
Function cgtase_cont code:
function dx=cgtase_model(t,x,ux,K1,H,Ea,R,Ed,dH,dS,dCp,Temp0,rol,N,Di,V,Vs,...
K2,Kd,Clx,qO2,YH,Ah,Cp,U,uxp,I,K,Tempj,pm)
dx = zeros(6,1);
%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 Ag
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)
%state variables
u = (ux*x(2))/((pm(1)+x(2))*(1+K1/H+H/K2)*(x(2)+pm(5)+(x(2)^2)/pm(1))*(pm(9)+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/pm(7)+pm(6))*(-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+pm(2))*u*x(1)-K*x(5));
dx(6) = (pm(8)*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);
The outcome by the ode45 using the guess value on parameter:
The fitness value graph by using GA solver to run the function ga_error
%

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

回答（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