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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 1stopt进行复杂动力学ode多参数拟合 求助@dingd
31/1返回列表
查看: 845  |  回复: 2
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

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
回复此楼

» 猜你喜欢

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

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

感谢参与，应助指数 +1
要么自己先写成1stOpt的代码形式，要么用文本、图片把原问题描述清楚。一堆Matlab代码，连看带猜费尽。
赞一下 回复此楼
2楼2014-09-25 09:01:27
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lrk7597

新虫 (初入文坛)
引用回帖:
1楼: Originally posted by winterhao at 2014-09-24 16:08:55
您好 我想求助您帮我用高版本1stopt进行一组参数拟合 我在matlab里用lsqnonlin每次优化出来的参数都是我给的初值，并没有实现优化功能，一直是local minimum possible。我没用过1stopt，所以我现在只能把matlab cod ...

我只想说大哥 能不能把软件分享一下 1101461558@qq.com
赞一下 回复此楼
3楼2018-04-19 11:35:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 ICPP2025 的主题更新
31/1返回列表
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭