应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 遗传算法优化微分方程中的参数
51/1返回列表
查看: 2390  |  回复: 8
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

For_study

金虫 (小有名气)

木虫

[求助] 遗传算法优化微分方程中的参数 已有1人参与

%拟合微分方程参数
%直接拟合反应的活化能和指前因子
function mynewtry4
clear;
clc;

%遗传算法
lb1=5*ones(1,22); lb2=-100000*ones(1,22); lb=[lb1,lb2];  ub1=100*ones(1,22); ub2=-1000*ones(1,22); ub=[ub1,ub2];
options = gaoptimset('Generations',1000,'StallGenLimit',50,...
     'StallTimeLimit',50,'TolFun',1e-12,'TolCon',1e-12,'MutationFcn',@mutationadaptfeasible);
[h0,fval,exitflag,reason,output,final_pop]=ga(@my_funtest1,44,options);

%输出参数
fprintf('\n\n遗传算法的估计数值:\n');
disp(h0);


%构造适应度函数my_funtest
function yfit=my_funtest1(h)
    tspan1=0:0.05:0.5;
    y01=[0.561,0.298,0.141,0 0 0 0 0];
[t,ycal1]=ode45(@myfun1,tspan1,y01,[],h);

y02=ycal1(11,;
tspan2=0.5:0.05:1;
[t,ycal2]=ode45(@myfun2,tspan2,y02,[],h);
ycalculate=ycal2(11,
yreal2=[0.0587 0 0 0.1988 0.3889 0.2309 0.0397 0.083];
for i=4:8
   ff(i)=(ycalculate(i)-yreal2(i))^2;
end
yfit=sum(ff)
end

function dy = myfun1(t,y,h)

% 动力学微分方程
% 共分为8集总，原料（饱和分SS为1，芳香分SA为2，胶质沥青质SR为3），柴油（D）为4,汽油（G）为5，液化气（LPG）为6，干气（Gas）为7，焦炭（C）为8，
% 由于遗传算法只能返回向量，所以反应常数不能是矩阵，只能是向量，第一个数表示反应物，第二个数表示生成物
%反应常数k（1，4）=k（1），k（1，5）=k（2），k（1，6）=k（3），k（1，7）=k（4），k（1，8）=k（5），k（2，4）=k（6），
%  k（2，5）=k（7），k（2，6）=k（8），k（2，7）=k（9），k（2，8）=k（10），k（3，4）=k（11），k（3，5）=k（12），
%  k（3，6）=k（13），k（3，7）=k（14），k（3，8）=k（15），k（4，5）=k（16），k（4，6）=k（17），k（4，7）=k（18），
%  k（4，8）=k（19），k（5，6）=k（20），k（5，7）=k（21），k（5，8）=k（22）
%
%
A=1;N=1;
deact=1;
dens=13.357;  %单位是kg/m^3
Swh=90;   %单位s
Tem=788; R=8.314;

k=zeros(1,22);
k(1)=h(1)*exp(-h(23)/(R*Tem));  % R常数，8.314，Tem温度，h表示活化能23-44，J/mol   指前因子1-22   kg/(m^3.s)
k(2)=h(2)*exp(-h(24)/(R*Tem));
k(3)=h(3)*exp(-h(25)/(R*Tem));
k(4)=h(4)*exp(-h(26)/(R*Tem));
k(5)=h(5)*exp(-h(27)/(R*Tem));
k(6)=h(6)*exp(-h(28)/(R*Tem));
k(7)=h(7)*exp(-h(29)/(R*Tem));
k(8)=h(8)*exp(-h(30)/(R*Tem));
k(9)=h(9)*exp(-h(31)/(R*Tem));
k(10)=h(10)*exp(-h(32)/(R*Tem));
k(11)=h(11)*exp(-h(33)/(R*Tem));
k(12)=h(12)*exp(-h(34)/(R*Tem));
k(13)=h(13)*exp(-h(35)/(R*Tem));

k(14)=h(14)*exp(-h(36)/(R*Tem));
k(15)=h(15)*exp(-h(37)/(R*Tem));
k(16)=h(16)*exp(-h(38)/(R*Tem));
k(17)=h(17)*exp(-h(39)/(R*Tem));
k(18)=h(18)*exp(-h(40)/(R*Tem));
k(19)=h(19)*exp(-h(41)/(R*Tem));

k(20)=h(20)*exp(-h(42)/(R*Tem));
k(21)=h(21)*exp(-h(43)/(R*Tem));
k(22)=h(22)*exp(-h(44)/(R*Tem));

dy(1)=-(k(1)+k(2)+k(3)+k(4)+k(5))*y(1)*A*N*deact*dens/Swh;  % A重芳烃失活系数，N碱氮吸附失活系数，deact催化剂结焦失活系数，dens密度，Swh真实重时空速
dy(2)=-(k(6)+k(7)+k(8)+k(9)+k(10))*y(2)*A*N*deact*dens/Swh;
dy(3)=-(k(11)+k(12)+k(13)+k(14)+k(15))*y(3)*A*N*deact*dens/Swh;
dy(4)=(k(1)*y(1)+k(6)*y(2)+k(11)*y(3)-(k(16)+k(17)+k(18)+k(19))*y(4))*A*N*deact*dens/Swh;
dy(5)=(k(2)*y(1)+k(7)*y(2)+k(12)*y(3)+k(16)*y(4)-(k(20)+k(21)+k(22))*y(5))*A*N*deact*dens/Swh;
dy(6)=(k(3)*y(1)+k(8)*y(2)+k(13)*y(3)+k(17)*y(4)+k(20)*y(5))*A*N*deact*dens/Swh;
dy(7)=(k(4)*y(1)+k(9)*y(2)+k(14)*y(3)+k(18)*y(4)+k(21)*y(5))*A*N*deact*dens/Swh;
dy(8)=(k(5)*y(1)+k(10)*y(2)+k(15)*y(3)+k(19)*y(4)+k(22)*y(5))*A*N*deact*dens/Swh;
dy=dy';

end




function dy = myfun2(t,y,h)

% 动力学微分方程
% 共分为8集总，原料（饱和分SS为1，芳香分SA为2，胶质沥青质SR为3），柴油（D）为4,汽油（G）为5，液化气（LPG）为6，干气（Gas）为7，焦炭（C）为8，
% 由于遗传算法只能返回向量，所以反应常数不能是矩阵，只能是向量，第一个数表示反应物，第二个数表示生成物
%反应常数k（1，4）=k（1），k（1，5）=k（2），k（1，6）=k（3），k（1，7）=k（4），k（1，8）=k（5），k（2，4）=k（6），
%  k（2，5）=k（7），k（2，6）=k（8），k（2，7）=k（9），k（2，8）=k（10），k（3，4）=k（11），k（3，5）=k（12），
%  k（3，6）=k（13），k（3，7）=k（14），k（3，8）=k（15），k（4，5）=k（16），k（4，6）=k（17），k（4，7）=k（18），
%  k（4，8）=k（19），k（5，6）=k（20），k（5，7）=k（21），k（5，8）=k（22）
%
%
A=1;N=1;
deact=1;
dens=13.357;Swh=90;
Tem=761; R=8.314;

k=zeros(1,22);
k(1)=h(1)*exp(-h(23)/(R*Tem));  % R常数，8.314，Tem温度，h表示活化能1-22，之前因子23-44
k(2)=h(2)*exp(-h(24)/(R*Tem));
k(3)=h(3)*exp(-h(25)/(R*Tem));
k(4)=h(4)*exp(-h(26)/(R*Tem));
k(5)=h(5)*exp(-h(27)/(R*Tem));
k(6)=h(6)*exp(-h(28)/(R*Tem));
k(7)=h(7)*exp(-h(29)/(R*Tem));
k(8)=h(8)*exp(-h(30)/(R*Tem));
k(9)=h(9)*exp(-h(31)/(R*Tem));
k(10)=h(10)*exp(-h(32)/(R*Tem));
k(11)=h(11)*exp(-h(33)/(R*Tem));
k(12)=h(12)*exp(-h(34)/(R*Tem));
k(13)=h(13)*exp(-h(35)/(R*Tem));

k(14)=h(14)*exp(-h(36)/(R*Tem));
k(15)=h(15)*exp(-h(37)/(R*Tem));
k(16)=h(16)*exp(-h(38)/(R*Tem));
k(17)=h(17)*exp(-h(39)/(R*Tem));
k(18)=h(18)*exp(-h(40)/(R*Tem));
k(19)=h(19)*exp(-h(41)/(R*Tem));

k(20)=h(20)*exp(-h(42)/(R*Tem));
k(21)=h(21)*exp(-h(43)/(R*Tem));
k(22)=h(22)*exp(-h(44)/(R*Tem));

dy(1)=-(k(1)+k(2)+k(3)+k(4)+k(5))*y(1)*A*N*deact*dens/Swh;  % A重芳烃失活系数，N碱氮吸附失活系数，deact催化剂结焦失活系数，dens密度，Swh真实重时空速
dy(2)=-(k(6)+k(7)+k(8)+k(9)+k(10))*y(2)*A*N*deact*dens/Swh;
dy(3)=-(k(11)+k(12)+k(13)+k(14)+k(15))*y(3)*A*N*deact*dens/Swh;
dy(4)=(k(1)*y(1)+k(6)*y(2)+k(11)*y(3)-(k(16)+k(17)+k(18)+k(19))*y(4))*A*N*deact*dens/Swh;
dy(5)=(k(2)*y(1)+k(7)*y(2)+k(12)*y(3)+k(16)*y(4)-(k(20)+k(21)+k(22))*y(5))*A*N*deact*dens/Swh;
dy(6)=(k(3)*y(1)+k(8)*y(2)+k(13)*y(3)+k(17)*y(4)+k(20)*y(5))*A*N*deact*dens/Swh;
dy(7)=(k(4)*y(1)+k(9)*y(2)+k(14)*y(3)+k(18)*y(4)+k(21)*y(5))*A*N*deact*dens/Swh;
dy(8)=(k(5)*y(1)+k(10)*y(2)+k(15)*y(3)+k(19)*y(4)+k(22)*y(5))*A*N*deact*dens/Swh;
dy=dy';

end
end


问题叙述，遗传算法未设定lb和ub时，程序可以正常运行，获得结果，当设定lb和ub后，程序依然可以运行（不会报错），但是只能运算一次优化后的结果（设置断点时发现）。遗传算法好像没有继续优化下去，希望有大神可以帮忙解答一下问题。。。。谢谢。。。。
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
第一颗纽扣扣错了。。。。
1楼 2015-12-08 20:52:11
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

For_study

金虫 (小有名气)

木虫
引用回帖:
4楼: Originally posted by dingd at 2015-12-09 16:06:39
换1stOpt试试，微分方程拟合简单好用的多。

参数太多，估计只有付费版的才能解决问题。。。。
赞一下 回复此楼
第一颗纽扣扣错了。。。。
8楼2015-12-11 15:05:32
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 9 个回答

For_study

金虫 (小有名气)

木虫
没人吗。。。。。。。。。。？
回复此楼
第一颗纽扣扣错了。。。。
2楼2015-12-09 14:14:41
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

For_study

金虫 (小有名气)

木虫
自己再顶一下。。。。。
回复此楼
第一颗纽扣扣错了。。。。
3楼2015-12-09 14:15:00
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

感谢参与，应助指数 +1
换1stOpt试试，微分方程拟合简单好用的多。
赞一下 回复此楼
4楼2015-12-09 16:06:39
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 9 个回答
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 085602 289分求调剂 +5 WWW西西弗斯 2026-03-24 5/250 2026-03-24 18:51 by jhhcooi
[考研] 274求调剂 +5 顾九笙要谦虚 2026-03-24 5/250 2026-03-24 18:43 by jhhcooi
[考研] 0854 考研调剂 招生了！AI 方向 +5 pk3725069 2026-03-19 17/850 2026-03-24 17:30 by zhouxuan..
[考研] 0703化学调剂，求导师收 +6 天天好运来上岸� 2026-03-24 6/300 2026-03-24 17:19 by dick_runner
[考研] 299求调剂 +7 某某某某位 2026-03-21 7/350 2026-03-24 15:24 by cuifj
[考研] 277分求调剂，跨调材料 +3 考研调剂lxh 2026-03-24 3/150 2026-03-24 13:52 by JourneyLucky
[考研] 305分求调剂（食品工程） +5 Sxy112 2026-03-21 7/350 2026-03-24 12:27 by 544594351
[考研] 086003食品工程求调剂 +4 淼淼111 2026-03-24 4/200 2026-03-24 11:53 by 544594351
[考研] 求材料，环境专业调剂 +3 18567500178 2026-03-18 3/150 2026-03-23 23:50 by 热情沙漠
[考研] 327求调剂 +5 prayer13 2026-03-23 5/250 2026-03-23 22:11 by 星空星月
[考研] 293求调剂 +3 涛涛Wjt 2026-03-22 5/250 2026-03-22 22:21 by jiangpengfei
[考研] 315分，诚求调剂，材料与化工085600 +3 13756423260 2026-03-22 3/150 2026-03-22 20:11 by edmund7
[考研] 一志愿华中农业071010，总分320求调剂 +5 困困困困坤坤 2026-03-20 6/300 2026-03-22 17:41 by hxsm
[考研] 324求调剂 +6 lucky呀呀呀鸭 2026-03-20 6/300 2026-03-22 16:01 by ColorlessPI
[考研] 材料 271求调剂 +5 展信悦_ 2026-03-21 5/250 2026-03-21 17:29 by 学员8dgXkO
[考研] 求调剂 +3 .m.. 2026-03-21 4/200 2026-03-21 16:25 by barlinike
[考研] 265求调剂 +12 梁梁校校 2026-03-19 14/700 2026-03-21 13:38 by lature00
[考研] 一志愿武理材料305分求调剂 +6 想上岸的鲤鱼 2026-03-18 7/350 2026-03-21 01:03 by JourneyLucky
[考研] 一志愿中海洋材料工程专硕330分求调剂 +8 小材化本科 2026-03-18 8/400 2026-03-20 23:16 by JourneyLucky
[考研] 353求调剂 +3 拉钩不许变 2026-03-20 3/150 2026-03-20 19:56 by JourneyLucky
信息提示
关闭
请填处理意见
关闭