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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » 动力学参数拟合
41/1返回列表
查看: 870  |  回复: 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返回列表
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭