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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 遗传算法优化微分方程中的参数
91/1返回列表
查看: 2359  |  回复: 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

金虫 (小有名气)

木虫
没人吗。。。。。。。。。。？
回复此楼
第一颗纽扣扣错了。。。。
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的回帖

yufeishu

银虫 (小有名气)
然而并看不懂
回复此楼
5楼2015-12-10 10:00:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ybkooo

至尊木虫 (著名写手)

^^
好多笑脸.
回复此楼
abcd
6楼2015-12-10 21:36:38
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

For_study

金虫 (小有名气)

木虫
引用回帖:
6楼: Originally posted by ybkooo at 2015-12-10 21:36:38
好多笑脸.

笑脸是因为笑脸的代码和我打出来的一样，所以在网页上显示为笑脸，其实笑脸的地方是：）。。。。。。。
赞一下 回复此楼
第一颗纽扣扣错了。。。。
7楼2015-12-11 15:04:45
已阅   回复此楼   关注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的回帖

For_study

金虫 (小有名气)

木虫
自己再给自己加把油。。。。
赞一下 回复此楼
第一颗纽扣扣错了。。。。
9楼2015-12-11 15:06:01
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 For_study 的主题更新
91/1返回列表
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[博后之家] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +6 3dfhjxgsh7 2026-02-22 9/450 2026-02-23 07:49 by w4l55oybr1
[考博] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +4 khieu8v8m0 2026-02-22 4/200 2026-02-23 06:46 by jsjzfl
[论文投稿] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +3 w89i99eaeh 2026-02-22 4/200 2026-02-23 06:36 by w4l55oybr1
[公派出国] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +3 khieu8v8m0 2026-02-22 5/250 2026-02-23 06:29 by w4l55oybr1
[硕博家园] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +4 khieu8v8m0 2026-02-22 8/400 2026-02-23 06:24 by w4l55oybr1
[考研] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +4 usprnugpzw 2026-02-21 10/500 2026-02-23 04:58 by 5jlh3qtdvx
[论文投稿] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +3 khieu8v8m0 2026-02-22 6/300 2026-02-23 02:08 by 5jlh3qtdvx
[考博] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +5 3dfhjxgsh7 2026-02-22 6/300 2026-02-23 02:04 by 5jlh3qtdvx
[教师之家] 版面费该交吗 +7 苹果在哪里 2026-02-22 8/400 2026-02-22 22:37 by otani
[基金申请] 基金正文30页指的是报告正文还是整个申请书 +5 successhe 2026-02-16 6/300 2026-02-22 21:38 by 山西悬空寺空悬�
[基金申请] 面上可以超过30页吧？ +4 阿拉贡aragon 2026-02-22 4/200 2026-02-22 21:22 by 山西悬空寺空悬�
[教师之家] 为什么中国大学教授们水了那么多所谓的顶会顶刊，但还是做不出宇树机器人？ +5 欢乐颂叶蓁 2026-02-21 5/250 2026-02-22 21:15 by 山西悬空寺空悬�
[论文投稿] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +4 usprnugpzw 2026-02-21 6/300 2026-02-22 19:48 by w89i99eaeh
[考研] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +3 3dfhjxgsh7 2026-02-22 4/200 2026-02-22 16:52 by khieu8v8m0
[找工作] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +3 usprnugpzw 2026-02-22 3/150 2026-02-22 16:37 by khieu8v8m0
[公派出国] 售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急 +3 usprnugpzw 2026-02-21 4/200 2026-02-22 16:27 by khieu8v8m0
[基金申请] “人文社科而论，许多学术研究还没有达到民国时期的水平” +4 苏东坡二世 2026-02-18 5/250 2026-02-22 16:07 by liangep1573
[基金申请] 什么是人一生最重要的？ +4 瞬息宇宙 2026-02-21 4/200 2026-02-22 11:44 by huagongfeihu
[基金申请] 今年春晚有几个节目很不错，点赞！ +11 瞬息宇宙 2026-02-16 12/600 2026-02-21 21:14 by lq493392203
[基金申请] 体制内长辈说体制内绝大部分一辈子在底层，如同你们一样大部分普通教师忙且收入低 +9 瞬息宇宙 2026-02-20 12/600 2026-02-21 10:39 by 欢乐颂叶蓁
信息提示
关闭
请填处理意见
关闭