| 查看: 581 | 回复: 0 | |||
[交流]
参数辨识代码一直在运行，不出结果，
|
|
这个是程序， function KineticsEst1_Int %用积分法进行SCR氨泄漏和氨覆盖度反应速率分析进行参数拟合 clear all; clc k0=[27.8 0.514 15.2 0.07 8.53e-13]; %参数初值 lb=[0 0 0 0 0]; %参数下限 ub=[+inf +inf +inf +inf +inf]; %参数上限 x0=[0 0 298.15]; %氨泄漏,氨覆盖度,下游温度初始值 KineticsData; yexp = ExpData(:,2:3); % yexp: 实验数据(氨泄漏x1 下游温度x2) % 使用函数fmincon()进行参数估计 [k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp); fprintf('\n使用函数fmincon()估计得到的参数值为:\n') fprintf('\tk1 = %.4f\n',k(1)) fprintf('\tk2 = %.4f\n',k(2)) fprintf('\tk3 = %.4f\n',k(3)) fprintf('\tk4 = %.4f\n',k(4)) fprintf('\tk5 = %.4f\n',k(5)) fprintf(' The sum of the squares is: %.1e\n\n',fval) k_fmincon = k; % 使用函数lsqnonlin()进行参数估计 [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp); ci = nlparci(k,residual,jacobian); fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n') Output % 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计 k0 = k_fmincon; [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp); ci = nlparci(k,residual,jacobian); fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n') Output % ------------------------------------------------------------------ function f = ObjFunc4Fmincon(k,x0,yexp) tspan = [1 20 40 60 80 100 140 180 220 280 320 360 400 500 600]; [t x] = ode45(@KineticEqs,tspan,x0,[],k); y(:,1) = x(:,1); y(:,2) = x(:,2); f = sum((y(:,1)-yexp(:,1)).^2)+ sum((y(:,2)-yexp(:,2)).^2); % ------------------------------------------------------------------ function f = ObjFunc4LNL(k,x0,yexp) tspan = [1 20 40 60 80 100 140 180 220 280 320 360 400 500 600]; [t x] = ode45(@KineticEqs,tspan,x0,[],k); y(:,1) = x(:,1); y(:,2) = x(:,2); f1 = y(:,1) - yexp(:,1)+y(:,2) - yexp(:,2); f = [f1;f2]; % ------------------------------------------------------------------ function dxdt = KineticEqs(t,x,k) dxdt = ... [ ( 1/(0.81*0.0000386)*5.364/1000000-x(1)*((287/101325)*(1/(0.81*0.0000386)*(1.4/3600)*x(3)+7.3*k(1)*sqrt(x(2))*(1-x(3))))+7.3*k(2)*exp(-k(3)/x(2))*x(3) ) ( k(4)*(1.4/3600)*(473.15-x(2))-k(5)*((x(2))^4-(298.15)^4) ) ( (x(1)*k(1)*sqrt(x(2))*(1-x(3)))-k(2)*exp(-k(3)/x(2))*x(3) ) ]; 这个是数据： %SCRKineticsData % 动力学数据: % t_时间 x1_氨泄漏_mol/m3 x2_温度 ExpData = ... [ 1 0 298.15 20 1.78231E-16 305.1437153 40 2.89597E-15 391.847062 60 6.1406E-13 458.2637625 80 7.83266E-13 472.3382312 100 1.68638E-12 473.6522317 140 1.79014E-11 473.7325518 180 2.18536E-10 473.7330945 220 2.50998E-09 473.7331684 280 9.70709E-08 473.7331458 320 1.11857E-06 473.7328618 360 1.27677E-05 473.7321413 400 0.000144696 473.7250369 500 0.009942007 473.1759208 600 0.010292012 473.149977 ] 请大神指点一下  |
回复此楼» 本帖@通知
» 猜你喜欢
【征稿】Springer Nature热电专题 Discover Materials / Discover Chemistry
已经有47人回复
-
25届硕士申请，材料计算方向
已经有2人回复
-
物理学I论文润色/翻译怎么收费?
已经有273人回复
-
课题组招1-2名硕士
已经有2人回复
-
机时越来越贵之后，我把一部分 DFT 流程搬到了 AI 计算平台上
已经有0人回复
-
斯坦福大学与Elsevier世界前2%学者新排名查询
已经有1人回复
-
27申博推荐：南京航空航天大学国际前沿院光学方向招收博士生！
已经有19人回复
-
华南师范大学+双一流学科物理学+广东省杰青团队2027年博士招生+机器学习驱动材料设计
已经有1人回复
-
apl投稿状态
已经有2人回复
-
华南师范大学+双一流学科物理学+广东省杰青团队2027年博士招生+机器学习驱动材料设计
已经有1人回复
-
日内瓦国际发明展申报完整流程（科研人员通用指南）
已经有10人回复
