24小时热门版块排行榜

返回列表

1064126551

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 916
帖子: 15
在线: 25.4小时
虫号: 3099587
注册: 2014-03-29
专业: 分离过程

[求助] 动力学参数拟合

原料是：甲醇、CO、氧气。主反应：甲醇、CO、氧气合成DMC，副反应：CO和氧气合成CO2
任务：拟合出生成DMC和CO2的动力学参数（指前因子、活化能、级数）
动力学方程为dy1/dw=f(y1、y2)和dy2/dw=f(y1、y2),具体形式见附图。
整理后的实验数据：N0(原料初始的摩尔数)、P（压力）、t(温度)、x10(甲醇初始的摩尔分数)、x20(CO初始的摩尔分数)、x30(O2初始的摩尔分数)、y1(DMC产物的摩尔分数)、y2(O2产物的摩尔分数)、W=0.375(催化剂的质量)。(每次的实验条件不同，原料经过固定床上催化剂的质量0.375g，反应后得到y1和y2)
目标函数：（y1拟合-y1）的平方和+（y2拟合-y2）的平方和最小
下面是自己编的代码：出现的问题，觉得是：第一：wspan有问题，实验条件不同，但是积分限是相同的。我不知道如何解决了
                                    第二：自己函数调用上，ode45调用的函数KineticsEqs(w,C,k,yexp)，C，yexp怎么表示？
希望大神帮忙，改一下代码，感激不尽
我的代码如下：
function k1k2k3k4k5k6k7k8k9
format long
clear all
clc
wspan=[0 0.375];%每次实验条件的积分限，方程为dy/dw, w=0.375
y0=[0 0];%实验每次的初值，即yDMC=0,yco2=0
k0=[8e3 2e4 -0.8 1  1  4e5 3e4 0.52 0.63];%参考文献上给的值
lb=[0 0 -2 -2 -2 0 0 -2 -2];%下限
ub=[+inf +inf 3 3 3 +inf +inf 3 3];%上限
data=...
[
0.3091 0.1000  413.1500 0.6240 0.2500 0.1270 0.0150 0.0180
0.3260 0.1000  413.1500 0.5000 0.3340 0.1660 0.0160 0.0200
0.2796 0.1000  413.3500 0.5310 0.3490 0.1210 0.0150 0.0180
0.2856 0.1000  416.5500 0.5710 0.2870 0.1420 0.0170 0.0330
0.2801 0.1000  423.1500 0.5300 0.3480 0.1220 0.0200 0.0400
0.3260 0.1000  423.2500 0.5000 0.3340 0.1660 0.0210 0.0410
0.3091 0.1000  423.3500 0.6240 0.2500 0.1270 0.0210 0.0420
0.2853 0.1000  424.3500 0.5720 0.2860 0.1420 0.0210 0.0420
0.3190 0.1000  426.8500 0.5110 0.3450 0.1440 0.0210 0.0450
0.2853 0.1000  432.3500 0.5720 0.2860 0.1420 0.0250 0.0540
0.3962 0.1000  433.1500 0.5240 0.3520 0.1240 0.0230 0.0530
0.3190 0.1000  434.4500 0.5110 0.3450 0.1440 0.0250 0.0540
0.3085 0.2000  413.3500 0.6250 0.2500 0.1250 0.0180 0.0310
0.2793 0.2000  414.2500 0.5310 0.3490 0.1200 0.0180 0.0320
0.2804 0.2000  422.8500 0.5290 0.3480 0.1230 0.0250 0.0500
0.2856 0.2000  423.4500 0.5710 0.2860 0.1430 0.0260 0.0510
0.3083 0.2000  423.9500 0.6250 0.2500 0.1240 0.0260 0.0500
0.3085 0.2000  432.8500 0.6250 0.2500 0.1250 0.0300 0.0590
0.2788 0.2000  434.3500 0.5320 0.3500 0.1180 0.0300 0.0580
0.2368 0.3000  411.1500 0.5010 0.3340 0.1650 0.0230 0.0430
0.2600 0.3000  423.4500 0.5700 0.2870 0.1420 0.0280 0.0510
0.2597 0.3000  433.1500 0.5710 0.2880 0.1410 0.0360 0.0630
0.3207 0.3000  433.3500 0.5550 0.3340 0.1110 0.0340 0.0610
0.2344 0.3000  435.4500 0.5060 0.3370 0.1570 0.0380 0.0660
];
% N0       P       T       x10（1初始量）x20    x30    y1       y2 最后采用ode45拟合结果y(1)、y(2).目标函数为：(y(1)-y1)的平方和+(y(2)-y2)的平方和最小
yexp=data(:,1:6);
[k,resnorm,residual,exitflag,output,lambda,jacobian]=...
lsqnonlin(@ObjFunc,k0,lb,ub,wspan,y0,yexp);
ci=nlparci(k,redidual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.9f ± %.9f\n',k(1),ci(1,2)-k(1))
fprintf('\tk2 = %.9f ± %.9f\n',k(2),ci(2,2)-k(2))
fprintf('\tk3 = %.9f ± %.9f\n',k(3),ci(3,2)-k(3))
fprintf('\tk1 = %.9f ± %.9f\n',k(4),ci(4,2)-k(1))
fprintf('\tk2 = %.9f ± %.9f\n',k(5),ci(5,2)-k(2))
fprintf('\tk3 = %.9f ± %.9f\n',k(6),ci(6,2)-k(3))
fprintf('\tk1 = %.9f ± %.9f\n',k(7),ci(7,2)-k(1))
fprintf('\tk2 = %.9f ± %.9f\n',k(8),ci(8,2)-k(2))
  fprintf('\tk2 = %.9f ± %.9f\n',k(9),ci(9,2)-k(2))
fprintf('  The sum of the squares is: %.9e\n\n',resnorm)
function f =ObjFunc(k,tspan,x0,yexp)
      [w Xsim]=ode45(@KineticsEqs,tspan,x0,[],k);
      Xsim1=Xsim(:,1);
      Xsim2=Xsim(:,2);
      ysim(:,1)=Xsim1(2:end);
      ysim(:,2)=Xsim1(2:end);
      size(ysim(:,1));
      size(ysim(:,2));
      size(yexp(:,7));
      size(yexp(:,8));
f=[(ysim(:,1)-yexp(:,7));(ysim(:,2)-yexp(:,8))];
function  dCdw=KineticsEqs(w,C,k,yexp)
      N0=yexp(:,1);
      P=yexp(:,2);
      T=yexp(:,3);
      x(1)=yexp(:,4);
      x(2)=yexp(:,5);
      x(3)=yexp(:,6);
r1=k(1)*exp(-k(2)./T).*P.^(k(3)+k(4)+k(5))*...
(C(1)*(1.5*x(1)-2)+0.5*x(5)*x(1))^(k(3))*...
(C(1)*(1.5*x(2)-1)+x(5)*(0.5*x(2)-1)+x(2))^(k(4))*...
(C(1)*(1.5*x(3)-0.5)+x(5)*(0.3*x(3)-0.5)+x(3))^k(5);
r2=k(6)*exp(-k(7)./T).*P.^(k(8)+k(9))*...
(C(1)*(1.5*x(2)-1)+C(2)*(0.5*x(2)-1)+x(2))^k(8)*...
(C(1)*(1.5*x(3)-0.5)+C(2)*(0.3*x(3)-0.5)+x(3))^k(9);
dy1dw=(1+1.5*C(1)+0.5*C(2))*(r1*(1.5*C(1)+1)+0.5*C(2)*r2)./N0;
dy2dw=(1+1.5*C(1)+0.5*C(2))*(r2*(1+0.5*C(1))+1.5*r1*C(2))./N0;
dCdw=[dy1dw;dy2dw];

r1和r2的表达式.jpg