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

24小时热门版块排行榜    

Znn3bq.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » 动力学参数拟合
41/1返回列表
查看: 947  |  回复: 3
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

清辉12365

新虫 (初入文坛)

[求助] 动力学参数拟合

原料是：甲醇、CO、氧气。主反应：甲醇、CO、氧气合成DMC，副反应：CO和氧气合成CO2.
要干的活：拟合出生成DMC和CO2的动力学参数（指前因子、活化能、级数）
动力学方程为dy1/dw=f(y1、y2)和dy2/dw=f(y1、y2),具体形式见附图。
整理后的实验数据：N0(原料初始的摩尔数)、P（压力）、t(温度)、x1(甲醇初始的摩尔分数)、x2(CO初始的摩尔分数)、x3(O2初始的摩尔分数)、y1(DMC产物的摩尔分数)、y2(O2产物的摩尔分数)、W=0.375(催化剂的质量)。
关键难度在：主要是两个常微分方程、数据多不知如何进行运算。下面是自己编的程序，出现问题并且解决不上来，希望大神指教
function YYHJIA
clear all;clc
global t T P N0 X1 X2 X3 ;
t=[140 140 140.2 143.4 150 150.1 150.2 151.2 153.7 159.2 160 161.3 140.2 141.1 149.7 150.3 150.8 159.7 161.2 138 150.3 160 160.2 162.3];
%温度
T=t+273.15;  
P=[0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3];
%压力
N0=[0.3091 0.3260 0.2796 0.2856 0.2801 0.3260 0.3091 0.2853 0.3190 0.2853 0.3962 0.3190 0.3085 0.2793 0.2804 0.2856 0.3083 0.3085 0.2788 0.2368 0.26 0.2597 0.3207 0.2344];
%原料初始摩尔
x1=[0.624 0.5 0.5310 0.5710 0.53 0.5 0.624 0.572 0.511 0.572 0.524 0.511 0.625 0.531 0.529 0.571 0.625 0.625 0.532 0.501 0.57 0.571 0.555 0.506];
%原料甲醇占的摩尔分数
x2=[0.25 0.334 0.349 0.287 0.348 0.334 0.25 0.286 0.345 0.286 0.352 0.345 0.25 0.349 0.348 0.286 0.25 0.25 0.35 0.334 0.287 0.288 0.334 0.337];
%原料CO摩尔分数
x3=[0.127 0.166 0.121 0.142 0.122 0.166 0.127 0.142 0.144 0.142 0.124 0.144 0.125 0.12 0.123 0.143 0.124 0.125 0.118 0.165 0.142 0.141 0.111 0.157];
%原料O2的摩尔分数
Y1=[0.015 0.016 0.015 0.017 0.02 0.021 0.021 0.021 0.021 0.025 0.023 0.025 0.018 0.018 0.025 0.026 0.026 0.03 0.03 0.023 0.028 0.036 0.034 0.038];
%出口DMC的摩尔数
Y2=[0.018 0.02 0.018 0.033 0.04 0.041 0.042 0.042 0.045 0.054 0.053 0.054 0.031 0.032 0.05 0.051 0.05 0.059 0.058 0.043 0.051 0.063 0.061 0.066];
%出口CO的摩尔数
k0=[8e4 2e4 -0.5 0.5 1 4e5 3e4 0.5 0.5 ];%(指前因子、活化能、级数)初始量
lb=[0 0 -3 -3 -3 0  0 -3 -3 ];%下限
ub=[+inf +inf +inf +inf +inf +inf +inf +inf +inf];%上限
%使用lsqnonlin函数进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(@ObjFunc,k0,lb,ub,[],T,P,N0,X1,X2,X3,Y1,Y2);
ci=nlparci(k,residual,jacobian);
function f=ObjFunc(k,T,P,N0,X1,X2,X3,Y1,Y2)
wspan=[0 0.375];
Y0=[Y1' Y2']
[m,n]=size(Y0);
Ycal = zeros(m,n);
for i=1:m
[w yy]=ode45(@Equations,wspan,Y0(i,,[],k);
Ycal(i,=yy(end,;
end
f=Ycal-Y0;
function dydw=Equations(w,yy,k)
%如何表达呢？这么多数字
r1=k(1)*exp*(-k(2)./T).*P.^(k(3)+k(4)+k(5))*...
    (yy(1).*(1.5*x1-2)+0.5*yy(2).*x1).^(k(3))*...
    (yy(1).*(1.5*x2-1)+yy(2).*(0.5*x2-1)+x2).^(k(4))*...
    (yy(1).*(1.5*x3-0.5)+yy(2).*(0.3*x3-0.5)+x3).^k(5);
r2=k(6)*exp*(-k(7)./T).*P.^(k(8)+k(9))*...
    (yy(1).*(1.5*x2-1)+yy(2).*(0.5*x2-1)+x2).^k(8)*...
    (yy(1).*(1.5*x3-0.5)+yy(2).*(0.3*x3-0.5)+x3).^k(9);
dy1dw=(1+1.5*yy(1)+0.5*yy(2)).*(r1.*(1.5*yy(1)+1)+0.5*yy(1).*r2)./N0;
dy2dw=(1+1.5*yy(1)+0.5*yy(2)).*(r2.*(1+0.5*yy(2))+1.5*r1.*yy(2))./N0;
dydw=[dy1dw;dy2dw];
出现的问题如图：

动力学参数拟合
出现的问题.jpg


动力学参数拟合-1
方程形式.jpg


动力学参数拟合-2
方程具体形式2.jpg
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2015-03-04 16:31:29
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

信彼南山

木虫 (著名写手)
exp后面是不要*的
其余的太长，没看完
赞一下(1人) 回复此楼
2楼2015-03-05 01:28:52
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

清辉12365

新虫 (初入文坛)
引用回帖:
2楼: Originally posted by 信彼南山 at 2015-03-05 01:28:52
exp后面是不要*的
其余的太长，没看完

我也看到了，请问一下，还有什么错误吗？？麻烦了~~
赞一下 回复此楼
3楼2015-03-05 08:35:46
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

清辉12365

新虫 (初入文坛)
引用回帖:
2楼: Originally posted by 信彼南山 at 2015-03-05 01:28:52
exp后面是不要*的
其余的太长，没看完

function YYHJIA
clear all;clc
global T P N0 x1 x2 x3;
t=[140 140 140.2 143.4 150 150.1 150.2 151.2 153.7 159.2 160 161.3 140.2 141.1 149.7 150.3 150.8 159.7 161.2 138 150.3 160 160.2 162.3];
%温度
T=t+273.15;  
P=[0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3];
%压力
N0=[0.3091 0.3260 0.2796 0.2856 0.2801 0.3260 0.3091 0.2853 0.3190 0.2853 0.3962 0.3190 0.3085 0.2793 0.2804 0.2856 0.3083 0.3085 0.2788 0.2368 0.26 0.2597 0.3207 0.2344];
%原料初始摩尔
x1=[0.624 0.5 0.5310 0.5710 0.53 0.5 0.624 0.572 0.511 0.572 0.524 0.511 0.625 0.531 0.529 0.571 0.625 0.625 0.532 0.501 0.57 0.571 0.555 0.506];
%原料甲醇占的摩尔分数
x2=[0.25 0.334 0.349 0.287 0.348 0.334 0.25 0.286 0.345 0.286 0.352 0.345 0.25 0.349 0.348 0.286 0.25 0.25 0.35 0.334 0.287 0.288 0.334 0.337];
%原料CO摩尔分数
x3=[0.127 0.166 0.121 0.142 0.122 0.166 0.127 0.142 0.144 0.142 0.124 0.144 0.125 0.12 0.123 0.143 0.124 0.125 0.118 0.165 0.142 0.141 0.111 0.157];
%原料O2的摩尔分数
Y1=[0.015 0.016 0.015 0.017 0.02 0.021 0.021 0.021 0.021 0.025 0.023 0.025 0.018 0.018 0.025 0.026 0.026 0.03 0.03 0.023 0.028 0.036 0.034 0.038];
%出口DMC的摩尔数
Y2=[0.018 0.02 0.018 0.033 0.04 0.041 0.042 0.042 0.045 0.054 0.053 0.054 0.031 0.032 0.05 0.051 0.05 0.059 0.058 0.043 0.051 0.063 0.061 0.066];
%出口CO的摩尔数
k0=[8e4 2e4 -0.5 0.5 1 4e5 3e4 0.5 0.5 ];%(指前因子、活化能、级数)初始量
lb=[0 0 -3 -3 -3 0  0 -3 -3 ];%下限
ub=[+inf +inf +inf +inf +inf +inf +inf +inf +inf];%上限
%使用lsqnonlin函数进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(@ObjFunc,k0,lb,ub,[],T,P,N0,x1,x2,x3,Y1,Y2);
ci=nlparci(k,residual,jacobian);
function f=ObjFunc(k,T,P,N0,x1,x2,x3,Y1,Y2)
wspan=[0 0.375];
Y0=[Y1' Y2'];
X=[N0' P' T' x1' x2' x3' Y1' Y2']
[m,n]=size(Y0);
Ycal = zeros(m,n);
for i=1:m
[w yy]=ode45(@Equations,wspan,Y0(i,,[],k);
Ycal(i,=yy(end,;
end
f=Ycal-Y0;
function dydw=Equations(w,yy,k,N0,P,T,x1,x2,x3)
global T P N0 x1 x2 x3;
%如何表达呢？这么多数字
r1=k(1)*exp(-k(2)./T).*P.^(k(3)+k(4)+k(5))*...
    (yy(1).*(1.5*x1-2)+0.5*yy(2).*x1).^(k(3))*...
    (yy(1).*(1.5*x2-1)+yy(2).*(0.5*x2-1)+x2).^(k(4))*...
    (yy(1).*(1.5*x3-0.5)+yy(2).*(0.3*x3-0.5)+x3).^k(5);
r2=k(6)*exp(-k(7)./T).*P.^(k(8)+k(9))*...
    (yy(1).*(1.5*x2-1)+yy(2).*(0.5*x2-1)+x2).^k(8)*...
    (yy(1).*(1.5*x3-0.5)+yy(2).*(0.3*x3-0.5)+x3).^k(9);
dy1dw=(1+1.5*yy(1)+0.5*yy(2)).*(r1.*(1.5*yy(1)+1)+0.5*yy(1).*r2)./N0;
dy2dw=(1+1.5*yy(1)+0.5*yy(2)).*(r2.*(1+0.5*yy(2))+1.5*r1.*yy(2))./N0;
dydw=[dy1dw;dy2dw];
图像为出现的错误结果
动力学参数拟合-3
图1.jpg
赞一下 回复此楼
4楼2015-03-05 09:36:48
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 清辉12365 的主题更新
41/1返回列表
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 考研调剂 +18 硕星赴 2026-04-09 19/950 2026-04-10 02:08 by belle15
[考研] 297求调剂 +8 Kwgyz 2026-04-09 8/400 2026-04-09 23:22 by may_新宇
[考研] 083200 305分 求二轮调剂 不接受跨专业 +3 Claireyyyy 2026-04-09 3/150 2026-04-09 23:17 by xiaohe123
[考研] 一志愿211电子信息347求调剂 +3 554916 2026-04-03 4/200 2026-04-09 22:24 by 554916
[考研] 调剂 +19 2261744733 2026-04-08 19/950 2026-04-09 19:11 by vgtyfty
[考研] 316求调剂 +3 想读研究生( ?∵ 2026-04-07 3/150 2026-04-09 17:31 by wp06
[考研] 299求调剂 +6 ZVVZ13 2026-04-08 6/300 2026-04-09 15:48 by wp06
[考研] 一志愿211，化学310分，本科重点双非，求调剂 +13 努力奋斗112 2026-04-08 13/650 2026-04-08 21:17 by 学员tURuqU
[考研] 347材料专硕求调剂 +18 zj8215216 2026-04-06 18/900 2026-04-08 16:27 by luoyongfeng
[考研] 机械工程264学硕求调剂 +3 qiushangxian 2026-04-06 3/150 2026-04-08 01:53 by Linzejun
[考研] 304求调剂（085602，过四级，一志愿985） +25 化工人999 2026-04-04 26/1300 2026-04-07 22:06 by hemengdong
[考研] 22408 331分求调剂 +4 y__1 2026-04-06 4/200 2026-04-06 17:26 by 土木硕士招生
[考研] 求助 +3 卡卡东88 2026-04-06 4/200 2026-04-06 15:28 by going home
[考研] 332求调剂 +17 小小孟... 2026-04-05 18/900 2026-04-06 09:51 by 蓝云思雨
[考研] 一志愿北交大材料工程总分358求调剂 +6 cs0106 2026-04-05 6/300 2026-04-05 16:34 by imissbao
[考研] 085600调剂 +9 东照照照 2026-04-04 9/450 2026-04-05 13:44 by ujn_zhuj
[考研] 考研调剂 +11 小sun要好运 2026-04-04 11/550 2026-04-05 08:02 by qlm5820
[考研] 22408求调剂 354分 可跨专业 +3 hannnnnnn 2026-04-04 3/150 2026-04-04 14:35 by 土木硕士招生
[考研] 387求调剂 +4 爱吃片豆土 2026-04-03 5/250 2026-04-04 08:10 by 岸上的一条鱼
[考研] 数二英二348求调剂 +4 hxdzj1 2026-04-03 5/250 2026-04-03 21:25 by zhq0425
信息提示
关闭
请填处理意见
关闭