当前位置: 首页 > 计算模拟 >用matlab求动力学参数

用matlab求动力学参数

作者 自制奶酪
来源: 小木虫 150 3 举报帖子
+关注

我不知道MATLAB怎么输入求动力学函数的方程,抄了木虫大佬的内容,请大家帮我把它弄到能运行出来,动力学参数拟合问题,求助大家~~~



function odes_fit
format long
clear all
clc
k0=[0 0 0 0 0 0 0 0];%参数初值
lb=[0 0 0 0 0 0 0 0];ub=[+inf +inf +inf +inf +inf +inf +inf +inf];%lb、ub:参数下限和上限
x0=[0 0 0 0 0];

data=...
    [0.5         1             1.5            2            3           4           5           6
          18.7    6.5    2.7    2.1   1.7   1.5   1.2   0.2;
          48.9    31.3   15.3   6.9   2.2   1.8   1.7   0.5;
          19.3    34.7   41.2   40.9   37.2   31.7   26.1   19.4;
          9.7    19.6    28.7   34.5   38.4   38.2   36.3   33.4;
          0.2    2.6    5.7    9.3   17.1   24.8   30.5   36.7;
];
x0=data(1,2:end);
tspan = [data(:,1)'];
yexp = [data(2:end,2) data(2:end,3) data(2:end,4) data(2:end,5) data(2:end,6) data(2:end,7) data(2:end,8)];
%使用函数fmincon()进行参数估计
[k,fval,flag]=fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.9f \n',k(1))
fprintf('\tk2 = %.9f \n',k(2))
fprintf('\tk3 = %.9f \n',k(3))
fprintf('\tk4 = %.9f \n',k(4))
fprintf('\tk5 = %.9f \n',k(5))
fprintf('\tk6 = %.9f \n',k(6))
fprintf('The sum of the squares is:%.le\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
figure(1)

ts=0(max(tspan)-min(tspan))/100):max(tspan);
[ts, ys] = ode45(@KineticsEqs,ts,x0,[],k);
yy = [data(:,2) data(:,3) data(:,4) data(:,5) data(:,6)];
figure(1)
plot(ts,ys(:,1),'b',tspan,yy(:,1),'bo');
figure(2)
plot(ts,ys(:,2),'r',tspan,yy(:,2),'ro');
figure(3)
plot(ts,ys(:,3),'k',tspan,yy(:,3),'ko');
figure(4)
plot(ts,ys(:,4),'g',tspan,yy(:,4),'ko');
figure(5)
plot(ts,ys(:,4),'g',tspan,yy(:,4),'ko');
figure(6)
plot(ts,ys(:,6),'m',tspan,yy(:,6),'mo');
%legend('C1的计算值','C1的实验值','C2的计算值','C2的实验值','C3的计算值','C3的实验值','C4的计算值','C4的实验值','C5的计算值','C5的实验值','C6的计算值','Location','best');



function f = ObjFunc(k,tspan,x0,yexp)           % 目标函数
[t, Xsim] = ode45(@KineticsEqs,tspan,x0,[],k);
Xsim1=Xsim(:,1);
Xsim2=Xsim(:,2);
Xsim3=Xsim(:,3);
Xsim4=Xsim(:,4);
Xsim5=Xsim(:,5);
Xsim6=Xsim(:,6);

ysim(:,1) = Xsim1(2:end);
ysim(:,2) = Xsim2(2:end);
ysim(:,3) = Xsim3(2:end);
ysim(:,4) = Xsim4(2:end);
ysim(:,5) = Xsim5(2:end);
ysim(:,6) = Xsim6(2:end);



f = [(ysim(:,1)-yexp(:,1)) (ysim(:,2)-yexp(:,2)) (ysim(:,3)-yexp(:,3)) (ysim(:,5)-yexp(:,5)) ...
    (ysim(:,6)-yexp(:,6))];

function dCdt = KineticsEqs(t,C,k)              % ODE模型方程
C1=C(1);C2=C(2);C3=C(3);C4=C(4);C5=C(5);C6=C(6);
k1=k(1);k2=k(2);k3=k(3);k4=k(4);k5=k(5);k6=k(6);k7=k(7);
dC1dt = -k1*C1-k7*C1;
dC2dt = k1*C1-k2*C2-k5*C2;
dC3dt = k2*C2-k3*C3-k6*C3;
dC4dt = k3*C3-k4*C5;
dC5dt = k4*C5;
dC6dt = k5*C2+k6*C3+k7*C1;

dCdt = [dC1dt;dC2dt;dC3dt;dC4dt;dC5dt;dC6dt]; @beefly 返回小木虫查看更多

今日热帖
  • 精华评论
  • 独孤神宇

    这个代码不完整

    lsqnonlin 函数 ObjFunc4LNL 缺失

  • 自制奶酪

    引用回帖:
    2楼: Originally posted by 独孤神宇 at 2021-01-05 21:03:26
    这个代码不完整

    lsqnonlin 函数 ObjFunc4LNL 缺失

    您好,这个要怎么修改啊,求助

  • dingd

    参考下1stOpt的计算结果:
    Weighted Root of Mean Square Error (RMSE): 3.52215359203778
    Weighted Sum of Squared Residual: 434.194807406662
    Correlation Coef. (R): 0.967264353331916
    R-Square: 0.935600329226609
    Determination Coef. (DC): 0.698480780131114
    F-Statistic: -0.866225958773538

    Parameter                  Best Estimate
    --------------------        -------------
    k1        1.18206935637483
    k2        1.23804230026633
    k3        0.392064236557302
    k4        0.237399839174299
    k5        2.4090372987211E-18
    k6        1.44475738382174E-24
    k7        0.719127412563461
    f Initial Value         19.2237430015159,

猜你喜欢
下载小木虫APP
与700万科研达人随时交流
  • 二维码
  • IOS
  • 安卓