版块导航: 正在加载中...

登录注册

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

24小时热门版块排行榜

返回列表

feng198800

金虫 (正式写手)

应助: 15 (小学生)
金币: 860.2
帖子: 415
在线: 242.1小时
虫号: 700353

[交流] 动力学参数拟合

CO变换反应，CO+H2O=CO2+H2，用固定床管式积分反应器，要拟合动力学方程中的参数！动力学方程如下图。转化后的东西学方程为dy/dw=-22.4./v*k(1)*exp(-k(2)/R.*T).*yco^k(3).*(yh2oin-(ycoin-yco))^k(4).*(yco2in+(ycoin-yco))^k(5).*(yh2in+(ycoin-yco))^k(6).*(1-β)，其中β是逆反应程度，要拟合指前因子、反应活化能和四个反应级数（k1~k6）。整理后的数据是已知反应器进口一氧化碳、二氧化碳、氢气、水蒸气的湿基组成和反应器出口一氧化碳湿基组成，主要是一个微分方程dy/dw=f(y),其中f(y)中含有要拟合的6个参数，涉及到微分方程的初值问题和最小二乘法问题。自己编了一个m文件，可是运行不了，纠结好长时间也没有发现问题，想大家帮我看一下！代码如下
function kinetics
clear all
clc
global t T R kp ycoin yh2oin yco2in yh2in v
t=[448.69 449.03 450.18 450.71 450.47 430.01 430.31 430.03 430.17 427.29 410.13 410.08 410.37 410.00 409.96 390.00 390.00 389.99 390.00 390.00];%反应温度
T=t+273.15;%绝对温度K
R=8.314;%气体常数
kp=exp(4575./T-4.33);%平衡常数
ycoin=[0.09865 0.13296 0.18299 0.23883 0.13809 0.13918 0.15477 0.14327 0.16374 0.13077 0.10624 0.11253 0.19904 0.20900 0.19895 0.12235 0.17431 0.14712 0.18064 0.16864];%进口CO含量
yh2oin=[0.52349 0.45073 0.39759 0.33124 0.45811 0.30844 0.39573 0.51067 0.54215 0.32203 0.48650 0.54179 0.32176 0.42416 0.32176 0.41059 0.32107 0.50444 0.48270 0.31474];%进口H2O含量
yco2in=[0.09440 0.08103 0.07302 0.05747 0.09028 0.11332 0.11931 0.04328 0.05115 0.10654 0.06181 0.03400 0.13821 0.09236 0.13688 0.04807 0.07740 0.07265 0.09696 0.08260];%进口CO2含量
yh2in=[0.09363 0.14281 0.17864 0.22005 0.13013 0.20040 0.20556 0.10088 0.10997 0.20593 0.17479 0.13200 0.16644 0.09602 0.17151 0.14214 0.13702 0.18113 0.16127 0.13905];%进口H2含量
v=[21.3952 18.2224 15.0234 13.5257 19.5599 14.5256 18.8683 19.0081 20.6600 16.2300 18.4185 20.6734 13.9253 15.8444 13.9233 16.3679 13.9538 20.7244 17.0155 14.2337];%总原料流量
ycoout=[0.07304 0.08526 0.12500 0.16207 0.08445 0.10695 0.13233 0.10264 0.11417 0.09878 0.09192 0.08657 0.17122 0.17377 0.17204 0.10501 0.14913 0.13370 0.15282 0.14068];%出口CO含量
k0=[3.5e5 8.4e4 0.9 0.6 -0.5 -0.04]; % k(指前因子、活化能、反应级数)初始向量
lb=[1e5 3e4 0 0 -1 -0.1];
ub=[1e6 1.5e5 1 1 0.5 0.5];% k上下限

% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag]=lsqnonlin(@ObjFunc,k0,lb,ub,[],ycoin,ycoout);

% ------------------------------------------------------------------
function f=ObjFunc(k,ycoin,ycoout)    % 目标函数
wspan=[0 0.6091];       % 即催化剂装填量/g
ycow=zeros(1,length(ycoout));
[w yco] =ode45(@Euqations,wspan,ycoin,[],k);%4,5阶龙格库塔法求解

ycow=yco(length(yco));%出口CO含量

f=ycoout-ycow;

% --------此方程定义微分方程等号右端项
function dydw=Euqations(w,yco,k)       % w催化剂质量
global T R kp ycoin yh2oin yco2in yh2in v
dydw=-22.4./v*k(1)*exp(-k(2)/R.*T).*yco^k(3).*(yh2oin-(ycoin-yco))^k(4).*(yco2in+(ycoin-yco))^k(5).*(yh2in+(ycoin-yco))^k(6).*(1-(yco2in-(ycoin-yco).*yh2in-(ycoin-yco))./(kp.*yco.*yh2oin-(ycoin-yco)));% 反应速率

[ 来自科研家族化工工程工艺 ]

回复此楼

» 收录本帖的淘帖专辑推荐

matlab

» 本帖已获得的红花（最新10朵）

dingd

» 猜你喜欢

【复旦大学】二维材料方向招收2026年博士研究生1名已经有0人回复
北京纳米能源与系统研究所王中林院士/曹南颖研究员课题组2026级硕/博/博后招生已经有10人回复
物理学I论文润色/翻译怎么收费? 已经有211人回复
荷兰Utrecht University超快太赫兹光谱王海教授课题招收2026 CSC博士生已经有23人回复
反铁磁体中的磁性切换：两种不同的机制已成功可视化已经有0人回复
求标准粉末衍射卡号 ICDD 01-076-1802 已经有0人回复
新西兰Robinson研究所招收全奖PhD 已经有0人回复
石墨烯转移--二氧化硅衬底石墨烯已经有0人回复

» 本主题相关价值贴推荐，对您同样有帮助:

请教matlab反应动力学参数估计遇到的问题，谢谢已经有15人回复
matlab 求指点动力学方程拟合过程中导数的获取已经有10人回复
matlab拟合拟均相动力学参数已经有6人回复
最小二乘拟合初值运算时间已经有12人回复
吸附的动力学方程拟合已经有18人回复
固定床中气-固相脱水反应本征动力学已经有5人回复
【讨论】关于腐蚀速度与动力学参数已经有15人回复
origin，sigmaplot等软件自定义曲线拟合（主要是一些动力学方程的模拟）已经有2人回复
【求助】如何从极化曲线求动力学参数已经有7人回复

» 抢金币啦！回帖就可以得到:

查看全部散金贴

1楼 2012-06-01 08:58:25

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

oYo

木虫 (著名写手)

应助: 5 (幼儿园)
金币: 3040.8
帖子: 1589
在线: 42.6小时
虫号: 400888

★
feng198800(金币+2): 谢谢参与

不懂，帮顶，

回复此楼

3楼2012-06-01 09:38:45

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

csgt0

荣誉版主 (著名写手)

.*,.^不完全，行向量和列向量不分，请检查公式。

赞一下

回复此楼

8楼2012-06-01 10:18:18

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

计算强帖: 4
应助: 1641 (讲师)
金币: 15037.3
帖子: 3410
在线: 1223.5小时
虫号: 291104

★
feng198800(金币+2): 谢谢参与

送鲜花一朵

微分方程拟合这类问题试试1stOpt，应该是最简单好用的了。

赞一下

回复此楼

9楼2012-06-01 12:01:47

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

feng198800

金虫 (正式写手)

应助: 15 (小学生)
金币: 860.2
帖子: 415
在线: 242.1小时
虫号: 700353

引用回帖:

8楼: Originally posted by csgt0 at 2012-06-01 10:18:18
.*,.^不完全，行向量和列向量不分，请检查公式。

你好，.*,.^不完全——可是公式中有的不是矩阵相乘处，可以不用点乘除的。请问您说的行向量和列向量不分指的是？？谢谢了

赞一下

回复此楼

11楼2012-06-01 14:47:56

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

csgt0

荣誉版主 (著名写手)

★ ★
cenwanglai: 金币+2, 公式输入辛苦了~ 2012-06-02 09:45:35

引用回帖:

11楼: Originally posted by feng198800 at 2012-06-01 14:47:56
你好，.*,.^不完全——可是公式中有的不是矩阵相乘处，可以不用点乘除的。请问您说的行向量和列向量不分指的是？？谢谢了...

dydw=-22.4./v*k(1)*exp(-k(2)/R.*T).*yco^k(3).*(yh2oin-(ycoin-yco))^k(4).*(yco2in+(ycoin-yco))^k(5).*(yh2in+(ycoin-yco))^k(6).*(1-(yco2in-(ycoin-yco).*yh2in-(ycoin-yco))./(kp.*yco.*yh2oin-(ycoin-yco)));% 反应速率

比如：yco^k(3)得用点的吧
ycoin和yco就一个横一个竖

赞一下

回复此楼

12楼2012-06-01 15:25:45

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

zhdgzhdg

至尊木虫 (著名写手)

应助: 32 (小学生)
金币: 22701.4
帖子: 1750
在线: 973.9小时
虫号: 95566

★
feng198800(金币+2): 谢谢参与

改变参数初值，或将参数归一化，比如参数之间差3个数量级以上，拟合时，小的那个参数就不会被调节。

赞一下

回复此楼

14楼2012-06-03 09:45:44

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

zhdgzhdg

至尊木虫 (著名写手)

应助: 32 (小学生)
金币: 22701.4
帖子: 1750
在线: 973.9小时
虫号: 95566

建议先根据文献指定四个反应级数，得到指前因子和活化能的值做为初值，在进行6个参数的拟合。

赞一下

回复此楼

15楼2012-06-03 09:51:27

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

dbb627

荣誉版主 (著名写手)

★
feng198800(金币+2): 谢谢参与

方程写错了，针对每个T用ode算ycow，不要用向量

赞一下

回复此楼

18楼2012-06-04 16:20:19

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

feng198800

金虫 (正式写手)

应助: 15 (小学生)
金币: 860.2
帖子: 415
在线: 242.1小时
虫号: 700353

引用回帖:

18楼: Originally posted by dbb627 at 2012-06-04 16:20:19
方程写错了，针对每个T用ode算ycow，不要用向量

谢谢了，可是我的动力学方程是每次实验的初始条件不同，也就是ode求解时的初值不同（温度T,流量V,进口湿基组成都不同）呀？

赞一下

回复此楼

19楼2012-06-05 22:42:55

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

dbb627

荣誉版主 (著名写手)

引用回帖:

19楼: Originally posted by feng198800 at 2012-06-05 22:42:55
谢谢了，可是我的动力学方程是每次实验的初始条件不同，也就是ode求解时的初值不同（温度T,流量V,进口湿基组成都不同）呀？...

没关系，我指的是方程形式一样，在调用时分开，参考
http://muchong.com/bbs/viewthread.php?tid=4277851&page=1#pid5

CODE:

function f = ObjFunc(k,P0,Pi)           % 目标函数

[m,n] = size(P0);

Pcal = zeros(m,n);

tspan =[0 0.0066];         % 即Wc/F0,kg.h/mol

for i = 1:m

[t PP] = ode45(@Euqations,tspan,P0(i,:),[],k);

Pcal(i,:) = PP(end,:);

end

f= Pcal-Pi;

% ------------------------------------------------------------------

function dPdt = Euqations(t, P, k)        % here t = Wc / F0

global keq1 keq2

m=k(6); n=k(7);

% 反应速率

r1 = k(1)*(P(6)*P(7)*P(7) - P(1)/keq1);

r2 = k(1)*P(1)*P(1);

r3 = k(3)*P(1)^(4-m)*P(2)^(m-2);

r4 = k(4)*P(6)^n*P(8)^(2*n+1);

r5 = k(5)*(P(6)*P(8) - P(5)*P(7)/keq2);

dPAdt = r1-2*r2-(4-m)*r3;

dPBdt = r2-(m-2)*r3;

dPCdt = r3;

dPDdt = r4;

dPEdt = r5;

dPFdt = r1+n*r4-r5;

dPGdt = -2*r1-(2*n+1)*r4+r5;

dPHdt = r2+r3+n*r4-r5;

dPdt = [dPAdt; dPBdt;dPCdt;dPDdt;dPEdt;dPFdt;dPGdt;dPHdt];

[ Last edited by dbb627 on 2012-6-6 at 16:25 ]

赞一下

回复此楼

20楼2012-06-06 16:23:50

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

strive123123

禁虫 (文坛精英)

★
feng198800(金币+2): 谢谢参与

本帖内容被屏蔽

21楼2012-06-07 09:45:26

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

shu_huiqiang

铜虫 (小有名气)

应助: 1 (幼儿园)
金币: 197.6
帖子: 64
在线: 37.4小时
虫号: 1541710

★
feng198800(金币+2): 谢谢参与

帮顶爱学术的人

回复此楼

22楼2012-10-24 16:05:28

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

forcal

禁虫 (初入文坛)

★ ★ ★ ★
小木虫: 金币+0.5, 给个红包，谢谢回帖
dbb627: 金币+3, 谢谢应助！ 2012-11-12 09:04:32

本帖内容被屏蔽

23楼2012-11-11 16:51:37

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

forcal

禁虫 (初入文坛)

★ ★ ★
小木虫: 金币+0.5, 给个红包，谢谢回帖
dbb627: 金币+2, 谢谢应助！ 2012-11-15 20:04:45

本帖内容被屏蔽

24楼2012-11-12 17:00:36

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

dragonpu

木虫 (小有名气)

应助: 0 (幼儿园)
金币: 2303.5
帖子: 71
在线: 79.8小时
虫号: 2112311

★ ★ ★
小木虫: 金币+0.5, 给个红包，谢谢回帖
csgt0: 金币+2, 谢谢 2012-11-17 09:08:43

Optimization terminated: relative function value
changing by less than OPTIONS.TolFun.

k =

  1.0e+005 *

  Columns 1 through 4

1.000003573678879 0.762988537687594 0.000005843208004 0.000005981921222

  Columns 5 through 6

  -0.000000744621941  -0.000000996174141

resnorm =

5.657234679833130e-004

residual =

  -0.006805767672396
0.002960729069534
  -0.006244673627901
  -0.001713302178392
0.010027828622409
0.004669076015188
  -0.006971457219631
  -0.002080652235825
0.005813395104906
0.007827253485404
  -0.008382737359081
0.000669267176915
  -0.001261508186825
  -0.002006057821641
  -0.001995402177405
  -0.001123333955809
0.003628239382343
  -0.004182989575642
0.004948094412108
0.007850845335425

exitflag =

   3

赞一下

回复此楼

25楼2012-11-15 16:59:53

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖