| 查看: 941 | 回复: 2 | ||
winterhao铁杆木虫 (正式写手)
|
[求助]
1stopt进行复杂动力学ode多参数拟合 求助@dingd 已有1人参与
|
|
您好 我想求助您帮我用高版本1stopt进行一组参数拟合 我在matlab里用lsqnonlin每次优化出来的参数都是我给的初值，并没有实现优化功能，一直是local minimum possible。我没用过1stopt，所以我现在只能把matlab code贴出来，请您有时间帮我看一下。非常感谢 %拟合函数 function fitting_kn_kg_n_g_f0_f1_f2_f3 format long clear all; clc % initial parameter identification kn = 10^2.6; kg = 8.06e-2; n = 6.5 ; g = 6.5 ; f=1; k0 = [kn,kg,n,g,f]; %Concentration and PH %Experimental data % Time tspan = [0,600,900,1200,1500,1800,2100,2400,2700,3000]; % Concentration and PH data = [0,0.186,0.271,0.36,0.384,0.311,0.262,0.257,0.252,0.245;... 10.35,8.64,8.45,8.35,8.21,8.09,8,7.83,7.65,7.42]'; size(data); yexp = data; %initial values for KineticsEqs x0 = [0,0,3.3e-4,3.03e-11,0,0,1.65e-4,100,0,0,0,0,0]; % Boundary condition of fitting parameter lb = [0 0 0 0 0 ]; ub = [1 1 1 1 1 ]*1e3; % Call the optimizer: options = optimoptions(@lsqnonlin,'TolX',1e-12,'TolFun',1e-12,'MaxFunEvals',500); [k,resnorm,residual,exitflag,output,lambda,jacobian] ... = lsqnonlin(@ObjFunc,k0,lb,ub,options,tspan,x0,yexp); %ci = nlparci(k,residual,jacobian); %Optimized solution ts = [0,3000]; [ts ys] = ode15s(@KineticsEqs,ts,x0,[],k);% the simulation data base on ts [ts_p ys_p] = ode15s(@KineticsEqs,tspan,x0,[],k);% the simulation data base on tspan R2_Mg = 1-sum((yexp(:,1)-ys_p(:,7)-ys_p(:,9)).^2)./sum((yexp(:,1)-mean(ys_p(:,7)+ys_p(:,9))).^2); R2_PH = 1-sum((yexp(:,2)+log10(ys_p(:,4))).^2)./sum((yexp(:,2)-mean(-log10(ys_p(:,4)))).^2); fprintf('n\t R2_Mg = %.6f\n',R2_Mg); fprintf('n\t R2_PH = %.6f\n',R2_PH); figure(1) plot(tspan,yexp(:,1),'o',ts,ys(:,7)+ys(:,9)) %plot the data vs solution title('Parameter fitting of Mg2+ in the bulk,mol/L'); figure(2) PH=-log10(ys(:,4)); plot(tspan,yexp(:,2),'o',ts,PH) %plot the data vs solution title('Parameter fitting of PH'); fprintf('\n\nlsqnonlin():\n') fprintf('\tkn = %.16f \n',k(1)) fprintf('\tkg = %.16f \n',k(2)) fprintf('\tn = %.16f \n',k(3)) fprintf('\tg = %.16f \n',k(4)) fprintf('\tf0 = %.16f \n',k(5)) end %优化函数 %Call 'Objfun' function: function lsq = ObjFunc(k,tspan,x0,yexp) [t Xsim] = ode15s(@KineticsEqs,tspan,x0,[],k); %Concentration of Mg2+ Xsim1 = Xsim(:,7)+Xsim(:,9); %PH Xsim2 = -log10(Xsim(:,4)); ycal(:,1) = Xsim1; ycal(:,2) = Xsim2; size(ycal(:,1)); size(ycal(:,2)); lsq = [(ycal(:,1)-yexp(:,1)) (ycal(:,2)-yexp(:,2))]; end %动力学方程组（ODE） % the reaction of Calc hydroxide[Mg(OH)2] and Carbon dioxide function dxdt=KineticsEqs(t,x,k) dxdt = zeros(13,1); Qg = 1.67e-2; % L/s 1L/min N = 650; vl = 2.826; % vl is the dispersion volum L vg = 2.047*10^(-5)*N^(0.579)*Qg^0.63*vl; % vg is the gas volum in the gas-liquid mixture m3 %agitating speed rpm P = 1; % P is the atmospheric pressure atm pi = 3.14; kla = 0.0304; % kla is the mass transfer coefficient of CO2 s^(-1) %(9.6218e-5)*N^0.4316*Qg^0.3840*c^0.1117*E; kH = 0.035; % kh is the Henry's constant mol/(L*atm) M_CO2 = 44; % M_CO2 is the co2 molar mass g/mol M_MgOH = 58; % molar mass of Mg(OH)2 g/mol M_MgCO3 = 84; % molar mass of M_MgCO3 g/mol dd_CO2 = 1.8; % dd_CO2 is the density of CO2 dd_MgOH = 2.34e3; % dd_Mg(OH)2 is the density of Mg(OH)2 g/L dd_MgCO3 = 2.96e3; % dd_Mg(OH)2 is the density of MgCO3 g/L dd_mol_MgCO3 = dd_MgCO3/M_MgCO3; % mol/L ddslurry = 1e3; % density of slurry g/L % ddCO2 = 1960; ddCO2 is the density of CO2 c = 1; % volume fraction of CO2 in the inlet gas CO2in = dd_CO2/M_CO2*c;% CO2in is the concentration of CO2 entering the reactor mol/L % CO2in = ddCO2/MCO2*0.15 =1800/44*0.05 =6.75 mol/m3 %Ksp1 = 1.8e-11; % Ksp[Mg(OH)2] mol2/L2 Ksp2 = 3.5e-8; % Ksp[MgCO3] mol2/L2 k11 = 8.4e3; % L/(mol*s) k11 = 8.4*10^3 L/mol*s k12 = 2e-4; % s^(-1) k12 = 2.0*10^(-4) s^(-1) k21 = 6e9; % L/(mol*s) k21 = 6*10^(9) L/mol*s k22 = 1.2e6; % s^(-1) k22 = 1.2*10^(6) s^(-1) k31 = 1.4e11; % L/(mol*s) k31 = 1.4*10^(11) L/mol*s k32 = 1.3e-3; % mol/(L*s) k32 = 1.3*10^(-3) mol/L*s k41 = 2.4e-2; % s^(-1) k41 = 2.4*10^(-2) s^(-1) k42 = 5.7e2; % L/(mol*s) k42 = 5.7*10^4 L/mol*s k52 = 9e-3; % s^(-1) k52 = 9e-3 s^(-1) k51 = 1.88e6; %k52/Ksp2; % L/(mol*s) D_OH = 5.27e-7; %diffusion coefficeint of OH- dm2/s D_H = 9.31e-7; %diffusion coefficeint of H+ dm2/s D_Mg =7.06e-8; %7.06e-8; %diffusion coefficeint of Mg2+ dm2/s zk_OH = -1; zk_H = 1; zk_Mg = 2; C_OH = 3.3e-4; C_H = 3.03e-11; C_Mg = 1.65e-4; % concentration on the interface i d_OH = C_OH-x(3); d_H = C_H-x(4); d_Mg = C_Mg-x(7); % Ci-Cb A_OH = (C_OH+x(3))/2; A_H = (C_H+x(4))/2; A_Mg = (C_Mg+x(7))/2; % (Ci+Cb)/2 Ntot1 = 8.17e10; %total number of Mg(OH)2 100g = 8.17e10 %Ntot1 = 100/(3.14*d^3*2340/6) UL = 1.01e-1; % dynamic viscosity 1.01e-3 Pas = kg*s/m = 1.01e-1 g*s/dm VL = 1.01e-4; % kinematic viscosity 1.01e-6 m2/s = 1.01e-4 dm2/s eT = 14028; % total energy dissipation w/kg = m2/s3*100 = dm2/s3 dp = (6*x(8)/(Ntot1*dd_MgOH*3.14))^(1/3); %diameter of solid dm a = (3.3*x(8)^(2/3)*(3.14*Ntot1)^(1/3)/dd_MgOH^(2/3)/vl);%specific area if a > 0 % mass transfer Re = eT^(1/3)*dp^(4/3)/VL; Sc_OH = UL/(ddslurry*D_OH); Sc_H = UL/(ddslurry*D_H); Sc_Mg = UL/(ddslurry*D_Mg); Sh_OH = (2^5.8+(0.61*(Re^0.58)*(Sc_OH^(1/3)))^5.8)^(1/5.8); Sh_H = (2^5.8+(0.61*(Re^0.58)*(Sc_H^(1/3)))^5.8)^(1/5.8); Sh_Mg = (2^5.8+(0.61*(Re^0.58)*(Sc_Mg^(1/3)))^5.8)^(1/5.8); Ks_OH = Sh_OH*D_OH/dp; Ks_H = Sh_H*D_H/dp; Ks_Mg = Sh_Mg*D_Mg/dp; K = (zk_OH*Ks_OH*d_OH+zk_H*Ks_H*d_H+zk_Mg*Ks_Mg*d_Mg)/... (zk_OH^2*Ks_OH*A_OH+zk_H^2*Ks_H*A_H+zk_Mg^2*Ks_Mg*A_Mg); NOH = Ks_OH*(d_OH-A_OH*zk_OH*K); NH = Ks_H*(d_H-A_H*zk_H*K); NMg = Ks_Mg*(d_Mg-A_Mg*zk_Mg*K); % reaction kinetics r1 = k11*x(1)*x(3)-k12*x(5); r2 = k21*x(5)*x(3)-k22*x(6); r3 = k31*x(3)*x(4)-k32; r4 = k41*x(1)-k42*x(5)*x(4); r5 = k51*x(7)*x(6)-k52; % Crystallization kinetics % Parameter fitting!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% k(1)-k(4) cMg=max(x(9),0); S = max(cMg-sqrt(Ksp2),0); B = k(1)*(S^k(3)); G = k(2)*(S^k(4)); %current = 2*NMg+NH-NOH %f=k(5)%%%%%%%%%%%%%%%%%%%%%%%%% k(5) dxdt(1) = kla*(kH*M_CO2*P*x(2)/dd_CO2-x(1))*(vl+vg)/vl-r1-r4; dxdt(2) = (Qg*(CO2in-x(2))-kla*(kH*M_CO2*P*x(2)/dd_CO2-x(1))*(vl+vg))/vg; dxdt(3) = -r1-r2-r3+NOH*a*k(5); dxdt(4) = -r3+r4+NH*a*k(5); dxdt(5) = r1-r2+r4; dxdt(6) = r2-r5; dxdt(7) = NMg*a*k(5)-r5; dxdt(8) = -NMg*a*k(5)*vl*M_MgOH; dxdt(9) = r5-pi/6*(1e-21)*dd_mol_MgCO3*B-G/2*pi*dd_mol_MgCO3*x(12); dxdt(10) = B; % m0 dxdt(11) = G*x(10);%m1 dxdt(12) = 2*G*x(11);%m2 dxdt(13) = 3*G*x(12);%m3 end end |
回复此楼» 猜你喜欢
［调剂信息］211智能人工感知方向国家青年特聘专家课题组招收调剂研究生
已经有0人回复
-
[调剂信息]211智能人工感知方向国家青年特聘专家课题组招收调剂研究生
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有287人回复
-
0702一志愿吉大B区求调剂 本科期间发表一篇Sci
已经有2人回复
-
070200求调剂，一志愿某211，288分
已经有18人回复
-
法国博士后职位
已经有0人回复
-
重庆交大26年硕士生招生拟调剂通知已出！欢迎加入机器视觉与3D光学成像课题组。
已经有0人回复
-
广州大学光电信息工程专业调剂，招收物理学专业学生
已经有1人回复
-
0702一志愿吉大B区求调剂 本科期间发表一篇Sci
已经有0人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 1stopt拟合中出现的问题
已经有8人回复
- 呼叫版主，在线紧急求助，关于matlab中微分方程组参数拟合得问题！
已经有12人回复
- 1stopt拟合出参数为0，为什么？
已经有6人回复
- 使用matlab最优化方法拟合获得多个动力学参数中的问题
已经有4人回复
- 求助1stopt拟合动力学参数
已经有4人回复
- 求大神帮忙拟合一个非线性方程，求出模型参数
已经有15人回复
- 用1stOpt进行复数曲线拟合时，结果出错。急求教！！！
已经有5人回复
- 求助用matlab拟合动力学方程
已经有9人回复
- 求大神帮忙用1stOpt 拟合复数数据
已经有13人回复
- 通过实验数据拟合，求解公式中的参数
已经有14人回复
- 求高版本1stopt拟合，
已经有11人回复
- 用1stopt拟合非线性方程结果与其他软件拟合结果差异大
已经有5人回复
- 如何拟合得到动力学模型参数？？
已经有3人回复
- 动力学方程参数估计方法
已经有14人回复
- 动力学参数拟合
已经有26人回复
- matlab 拟合反应动力学参数结果很差。大家帮忙看一下
已经有14人回复
- 【求助】用matlab最优化方法进行参数拟合
已经有17人回复
- 【求助】拟合动力学方程求助
已经有13人回复
dingd
铁杆木虫 (职业作家)
- 计算强帖: 4
- 应助: 1641 (讲师)
- 金币: 15037.3
- 散金: 101
- 红花: 234
- 帖子: 3410
- 在线: 1223.7小时
- 虫号: 291104
- 注册: 2006-10-28
2楼2014-09-25 09:01:27
3楼2018-04-19 11:35:17
