24小时热门版块排行榜

返回列表

小小莓邮

铁虫 (小有名气)

应助: 0 (幼儿园)
金币: 0.2
帖子: 52
在线: 48.8小时
虫号: 1245029
注册: 2011-03-25
性别: GG
专业: 化学反应工程

[求助] 求助遗传算法对反应速率常数初步估算的问题

各位大神小弟先在此谢过
我要用matlab做反应速率常数的拟合，打算用遗传算法对常数先进性拟合再用最小二乘法拟合，但用遗传算法是程序进行不下去，请各位帮忙看一下，问题出在哪，急求。
数据：
kineticsA4data
% 动力学数据

mmol/m2)
% t(h) C    A    P    P
ExpData = ...
[ 0.00 25.2751 0    0    0
   0.25 22.9500 0.3201 1.0294 0.3103
   0.50 21.6138 0.4853 1.4568 0.6454
   0.75 21.2019 0.5185 1.8947 0.9272
   1.00 20.3590 0.6083 2.1379 1.1658
   1.25 20.0490 0.6995 2.3849 1.2432
   1.50 19.3850 0.7270 2.4681 1.5247
   1.75 19.1517 0.7706 2.6907 1.6407
   2.00 18.9174 0.8237 3.1265 1.7000
   2.25 18.1695 0.9126 3.3860 1.9000
   2.50 18.3046 0.9659 3.3133 2.0000
   2.75 16.7374 0.9784 4.1797 2.1000
   3.00 17.2217 1.0000 3.9054 2.3000];
c程序：
function ga_test_fuzamodel_feifor
clear all
clc
lb = [0  0  0  0  0 0 0 0 0];                      % 参数下限
ub = [+inf  +inf  +inf  +inf  +inf +inf +inf +inf +inf];
x0=[25.2751 0 0 0 0.0001 0.0001 0.0001];
kineticsA4data;
yexp=ExpData(:,2:5);
options = gaoptimset('Generations',1000,'StallGenLimit',300,...
'StallTimeLimit',50,'TolFun',1e-12,'TolCon',1e-12);
[k1,fva,reason,output,final_pop]=ga(@objfun,9,options);
options = gaoptimset('InitialPopulation',final_pop,'Generations',1000,'StallGenLimit',300,...
'StallTimeLimit',50,'TolFun',1e-12,'TolCon',1e-12);
[k2,fva,reason,output,final_pop2]=ga(@objfun,9,options);
fprintf('\n\n遗传算法的初始估计数值:\n');
fprintf('\n\t参数 a0 = %.9f',k2(1));
fprintf('\n\t参数 b0 = %.9f',k2(2));
fprintf('\n\t参数 c0 = %.9f',k2(3));
fprintf('\n\t参数 d0 = %.9f',k2(4));
fprintf('\n\t参数 d0 = %.9f',k2(5));
fprintf('\n\t参数 d0 = %.9f',k2(6));
fprintf('\n\t参数 d0 = %.9f',k2(7));
fprintf('\n\t参数 d0 = %.9f',k2(8));
fprintf('\n\t参数 d0 = %.9f',k2(9));
function fun=objfun(k)
tspan=[0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5 2.75 3];
x0=[25.2751 0 0 0 0.0001 0.0001 0.0001];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);
kineticsA4data;
yexp=ExpData(:,2:5);
ycal(:,1) = x(:,1);
ycal(:,2) = x(:,2);
ycal(:,3) = x(:,3);
ycal(:,4) = x(:,4)
fun = sum((ycal(:,1)-yexp(:,1)).^2) + sum((ycal(:,2)-yexp(:,2)).^2) ...
+ sum((ycal(:,3)-yexp(:,3)).^2) + sum((ycal(:,4)-yexp(:,4)).^2)
% ------------------------------------------------------------------
function dxdt = KineticEqs(t,x,k)
dxdt=[(-k(3)*x(1)*x(7)-k(4)*x(1)*x(6))
(k(7)*x(6))
(k(4)*x(1)*x(6)+k(5)*x(4)*x(6))
(k(3)*x(1)*x(7)-k(1)*x(4)-k(5)*x(4)*x(6))
(-k(2)*x(5)+k(3)*x(1)*x(7)+k(4)*x(1)*x(6)-2*k(9)*x(5)^2)
(k(1)*x(4)-k(4)*x(1)*x(6)-k(5)*x(4)*x(6)-k(7)*x(6)+2*k(6)*x(6)^2)
(k(2)*x(5)-k(3)*x(1)*x(7)+k(5)*x(4)*x(6)-2*(k(6)+k(8))*x(7)^2)];

回复此楼