常微分方程拟合matlab,求助matlab常微分方程参数拟合程序修改 - 计算模拟 - 小木虫 - 学术 科研 互动社区...

版主及各位高人，帮小弟看看这个程序吧

想拟合常微分方程的14个参数，从k1到k14.

别看表达式复杂，但是都有一样的重复结构单元

方程组和数据如下：

dy1/dt=-1/sw*(k1+k2+k3+k4+k5)*y1/(y1/644+y2/230+y3/115+y4/52+y5/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-y1-y2-y3-y4-y5)/6)^-0.385;

dy2/dt=1/sw*(k1*y1-(k8+k9+k10+k11)*y2)/(y1/644+y2/230+y3/115+y4/52+y5/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-y1-y2-y3-y4-y5)/6)^-0.385;

dy3/dt=1/sw*(k2*y1+k8*y2-(k12+k13+k14)*y3)/(y1/644+y2/230+y3/115+y4/52+y5/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-y1-y2-y3-y4-y5)/6)^-0.385;

dy4/dt=1/sw*(k3*y1+k9*y2+k12*y3)/(y1/644+y2/230+y3/115+y4/52+y5/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-y1-y2-y3-y4-y5)/6)^-0.385;

dy5/dt=1/sw*(k4*y1+k10*y2+k13*y3)/(y1/644+y2/230+y3/115+y4/52+y5/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-y1-y2-y3-y4-y5)/6)^-0.385;

sw=5,

t   y1   y2   y3   y4   y5

0  93.7   6.3  0    0    0

1  5.14  12.76  45.41  21.52  3.31

sw=10,

t   y1   y2   y3   y4   y5

0  93.7   6.3  0    0    0

1  5.36  13.20  46.94  20.90  3.17

sw=15,

t   y1   y2   y3   y4   y5

0  93.7   6.3  0    0    0

1  5.51  13.47  48.68  20.01  3.08

sw=20,

t   y1   y2   y3   y4   y5

0  93.7   6.3  0    0    0

1  5.78  14.03  48.16  19.34  3.04

sw=25,

t   y1   y2   y3   y4   y5

0  93.7   6.3  0    0    0

1  6.07  14.33  47.80  18.65  2.94

程序如下，但拟合的参数结果不理想，太大

function k1_k14

format long

clear all

clc

tspan = [0 1];

x0=[93.7 6.3 0 0 0 93.7 6.3 0 0 0 93.7 6.3 0 0 0 93.7 6.3 0 0 0 93.7 6.3 0 0 0];

k0 = [4 20 8 1.5 3 0.06 4 2 0.5 0.5 0.5 2 0.5 0.5];

lb = [0 0 0 0 0 0 0 0 0 0 0 0 0 0];

ub = [1 1 1 1 1 1 1 1 1 1 1 1 1 1]*1e5;

data=[5.14 12.76 45.41 21.52 3.31 5.36 13.20 46.94 20.90 3.17 5.51 13.47 48.68 20.01 3.08 5.78 14.03 48.16 19.34 3.04 6.07 14.33 47.80 18.65 2.94];

size(data);

yexp = data;

[k,fval,flag] =fmincon(@ObjFuncFmincon,k0,[],[],[],[],lb,ub,[],[],tspan,x0,yexp);

fprintf('\n\n使用函数fmincon()估计得到的参数值为:\n')

fprintf('\tk1 = %.16f \n',k(1))

fprintf('\tk2 = %.16f \n',k(2))

fprintf('\tk3 = %.16f \n',k(3))

fprintf('\tk4 = %.16f \n',k(4))

fprintf('\tk5 = %.16f \n',k(5))

fprintf('\tk6 = %16f \n',k(6))

fprintf('\tk7 = %.16f \n',k(7))

fprintf('\tk8 = %.16f \n',k(8))

fprintf('\tk9 = %.16f \n',k(9))

fprintf('\tk10 = %.16f \n',k(10))

fprintf('\tk11 = %.16f \n',k(11))

fprintf('\tk12 = %.16f \n',k(12))

fprintf('\tk13 = %.16f \n',k(13))

fprintf('\tk14 = %.16f \n',k(14))

fprintf('The sum of squares is: %.16f \n\n',fval)

[k,resnorm,residual,exitflag,output,lambda,jacobian] =lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);

ci = nlparci(k,residual,jacobian);

fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')

fprintf('\tk1 = %.16f \n',k(1))

fprintf('\tk2 = %.16f \n',k(2))

fprintf('\tk3 = %.16f \n',k(3))

fprintf('\tk4 = %.16f \n',k(4))

fprintf('\tk5 = %.16f \n',k(5))

fprintf('\tk6 = %16f \n',k(6))

fprintf('\tk7 = %.16f \n',k(7))

fprintf('\tk8 = %.16f \n',k(8))

fprintf('\tk9 = %.16f \n',k(9))

fprintf('\tk10 = %.16f \n',k(10))

fprintf('\tk11 = %.16f \n',k(11))

fprintf('\tk12 = %.16f \n',k(12))

fprintf('\tk13 = %.16f \n',k(13))

fprintf('\tk14 = %.16f \n',k(14))

fprintf('The sum of squares is: %.16f \n\n',resnorm)

function f = ObjFuncFmincon(k,tspan,x0,yexp) % fmincon目标函数

[t Xsim] = ode45(@KineticsEqs,tspan,x0,[],k);

Xsim=Xsim';

Xsim1=Xsim(:,1);

Xsim2=Xsim(:,2);

Xsim3=Xsim(:,3);

Xsim4=Xsim(:,4);

Xsim5=Xsim(:,5);

Xsim6=Xsim(:,6);

Xsim7=Xsim(:,7);

Xsim8=Xsim(:,8);

Xsim9=Xsim(:,9);

Xsim10=Xsim(:,10);

Xsim11=Xsim(:,11);

Xsim12=Xsim(:,12);

Xsim13=Xsim(:,13);

Xsim14=Xsim(:,14);

Xsim15=Xsim(:,15);

Xsim16=Xsim(:,16);

Xsim17=Xsim(:,17);

Xsim18=Xsim(:,18);

Xsim19=Xsim(:,19);

Xsim20=Xsim(:,20);

Xsim21=Xsim(:,21);

Xsim22=Xsim(:,22);

Xsim23=Xsim(:,23);

Xsim24=Xsim(:,24);

Xsim25=Xsim(:,25);

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);

ysim(:,7) = Xsim7(2:end);

ysim(:,8) = Xsim8(2:end);

ysim(:,9) = Xsim9(2:end);

ysim(:,10) = Xsim10(2:end);

ysim(:,11) = Xsim11(2:end);

ysim(:,12) = Xsim12(2:end);

ysim(:,13) = Xsim13(2:end);

ysim(:,14) = Xsim14(2:end);

ysim(:,15) = Xsim15(2:end);

ysim(:,16) = Xsim16(2:end);

ysim(:,17) = Xsim17(2:end);

ysim(:,18) = Xsim18(2:end);

ysim(:,19) = Xsim19(2:end);

ysim(:,20) = Xsim20(2:end);

ysim(:,21) = Xsim21(2:end);

ysim(:,22) = Xsim22(2:end);

ysim(:,23) = Xsim23(2:end);

ysim(:,24) = Xsim24(2:end);

ysim(:,25) = Xsim25(2:end);

size(ysim(:,1));

size(ysim(:,2));

size(ysim(:,3));

size(ysim(:,4));

size(ysim(:,5));

size(ysim(:,6));

size(ysim(:,7));

size(ysim(:,8));

size(ysim(:,9));

size(ysim(:,10));

size(ysim(:,11));

size(ysim(:,12));

size(ysim(:,13));

size(ysim(:,14));

size(ysim(:,15));

size(ysim(:,16));

size(ysim(:,17));

size(ysim(:,18));

size(ysim(:,19));

size(ysim(:,20));

size(ysim(:,21));

size(ysim(:,22));

size(ysim(:,23));

size(ysim(:,24));

size(ysim(:,25));

size(yexp(:,1));

size(yexp(:,2));

size(yexp(:,3));

size(yexp(:,4));

size(yexp(:,5));

size(yexp(:,6));

size(yexp(:,7));

size(yexp(:,8));

size(yexp(:,9));

size(yexp(:,10));

size(yexp(:,11));

size(yexp(:,12));

size(yexp(:,13));

size(yexp(:,14));

size(yexp(:,15));

size(yexp(:,16));

size(yexp(:,17));

size(yexp(:,18));

size(yexp(:,19));

size(yexp(:,20));

size(yexp(:,21));

size(yexp(:,22));

size(yexp(:,23));

size(yexp(:,24));

size(yexp(:,25));

f = sum((ysim(:,1)-yexp(:,1)).^2)+sum((ysim(:,2)-yexp(:,2)).^2)+sum((ysim(:,3)-yexp(:,3)).^2)+sum((ysim(:,4)-yexp(:,4)).^2)+sum((ysim(:,5)-yexp(:,5)).^2)+sum((ysim(:,6)-yexp(:,6)).^2)+sum((ysim(:,7)-yexp(:,7)).^2)+sum((ysim(:,8)-yexp(:,8)).^2)+sum((ysim(:,9)-yexp(:,9)).^2)+sum((ysim(:,10)-yexp(:,10)).^2)+sum((ysim(:,11)-yexp(:,11)).^2)+sum((ysim(:,12)-yexp(:,12)).^2)+sum((ysim(:,13)-yexp(:,13)).^2)+sum((ysim(:,14)-yexp(:,14)).^2)+sum((ysim(:,15)-yexp(:,15)).^2)+sum((ysim(:,16)-yexp(:,16)).^2)+sum((ysim(:,17)-yexp(:,17)).^2)+sum((ysim(:,18)-yexp(:,18)).^2)+sum((ysim(:,19)-yexp(:,19)).^2)+sum((ysim(:,20)-yexp(:,20)).^2)+sum((ysim(:,21)-yexp(:,21)).^2)+sum((ysim(:,22)-yexp(:,22)).^2)+sum((ysim(:,23)-yexp(:,23)).^2)+sum((ysim(:,24)-yexp(:,24)).^2)+sum((ysim(:,25)-yexp(:,25)).^2);

function f = ObjFunc(k,tspan,x0,yexp) % lsqnonlin目标函数

[t Xsim] = ode45(@KineticsEqs,tspan,x0,[],k);

Xsim=Xsim';

Xsim1=Xsim(:,1);

Xsim2=Xsim(:,2);

Xsim3=Xsim(:,3);

Xsim4=Xsim(:,4);

Xsim5=Xsim(:,5);

Xsim6=Xsim(:,6);

Xsim7=Xsim(:,7);

Xsim8=Xsim(:,8);

Xsim9=Xsim(:,9);

Xsim10=Xsim(:,10);

Xsim11=Xsim(:,11);

Xsim12=Xsim(:,12);

Xsim13=Xsim(:,13);

Xsim14=Xsim(:,14);

Xsim15=Xsim(:,15);

Xsim16=Xsim(:,16);

Xsim17=Xsim(:,17);

Xsim18=Xsim(:,18);

Xsim19=Xsim(:,19);

Xsim20=Xsim(:,20);

Xsim21=Xsim(:,21);

Xsim22=Xsim(:,22);

Xsim23=Xsim(:,23);

Xsim24=Xsim(:,24);

Xsim25=Xsim(:,25);

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);

ysim(:,7) = Xsim7(2:end);

ysim(:,8) = Xsim8(2:end);

ysim(:,9) = Xsim9(2:end);

ysim(:,10) = Xsim10(2:end);

ysim(:,11) = Xsim11(2:end);

ysim(:,12) = Xsim12(2:end);

ysim(:,13) = Xsim13(2:end);

ysim(:,14) = Xsim14(2:end);

ysim(:,15) = Xsim15(2:end);

ysim(:,16) = Xsim16(2:end);

ysim(:,17) = Xsim17(2:end);

ysim(:,18) = Xsim18(2:end);

ysim(:,19) = Xsim19(2:end);

ysim(:,20) = Xsim20(2:end);

ysim(:,21) = Xsim21(2:end);

ysim(:,22) = Xsim22(2:end);

ysim(:,23) = Xsim23(2:end);

ysim(:,24) = Xsim24(2:end);

ysim(:,25) = Xsim25(2:end);

size(ysim(:,1));

size(ysim(:,2));

size(ysim(:,3));

size(ysim(:,4));

size(ysim(:,5));

size(ysim(:,6));

size(ysim(:,7));

size(ysim(:,8));

size(ysim(:,9));

size(ysim(:,10));

size(ysim(:,11));

size(ysim(:,12));

size(ysim(:,13));

size(ysim(:,14));

size(ysim(:,15));

size(ysim(:,16));

size(ysim(:,17));

size(ysim(:,18));

size(ysim(:,19));

size(ysim(:,20));

size(ysim(:,21));

size(ysim(:,22));

size(ysim(:,23));

size(ysim(:,24));

size(ysim(:,25));

size(yexp(:,1));

size(yexp(:,2));

size(yexp(:,3));

size(yexp(:,4));

size(yexp(:,5));

size(yexp(:,6));

size(yexp(:,7));

size(yexp(:,8));

size(yexp(:,9));

size(yexp(:,10));

size(yexp(:,11));

size(yexp(:,12));

size(yexp(:,13));

size(yexp(:,14));

size(yexp(:,15));

size(yexp(:,16));

size(yexp(:,17));

size(yexp(:,18));

size(yexp(:,19));

size(yexp(:,20));

size(yexp(:,21));

size(yexp(:,22));

size(yexp(:,23));

size(yexp(:,24));

size(yexp(:,25));

f = [(ysim(:,1)-yexp(:,1)) (ysim(:,2)-yexp(:,2)) (ysim(:,3)-yexp(:,3)) (ysim(:,4)-yexp(:,4)) (ysim(:,5)-yexp(:,5)) (ysim(:,6)-yexp(:,6)) (ysim(:,7)-yexp(:,7)) (ysim(:,8)-yexp(:,8)) (ysim(:,9)-yexp(:,9)) (ysim(:,10)-yexp(:,10)) (ysim(:,11)-yexp(:,11)) (ysim(:,12)-yexp(:,12)) (ysim(:,13)-yexp(:,13)) (ysim(:,14)-yexp(:,14)) (ysim(:,15)-yexp(:,15)) (ysim(:,16)-yexp(:,16)) (ysim(:,17)-yexp(:,17)) (ysim(:,18)-yexp(:,18)) (ysim(:,19)-yexp(:,19)) (ysim(:,20)-yexp(:,20)) (ysim(:,21)-yexp(:,21)) (ysim(:,22)-yexp(:,22)) (ysim(:,23)-yexp(:,23)) (ysim(:,24)-yexp(:,24)) (ysim(:,25)-yexp(:,25))];

function dYdt = KineticsEqs(t,Y,k)

% ODE模型方程

Y1=Y(1);Y2=Y(2);Y3=Y(3);Y4=Y(4);Y5=Y(5);Y6=Y(6);Y7=Y(7);Y8=Y(8);Y9=Y(9);Y10=Y(10);Y11=Y(11);Y12=Y(12);Y13=Y(13);Y14=Y(14);Y15=Y(15);Y16=Y(16);Y17=Y(17);Y18=Y(18);Y19=Y(19);Y20=Y(20);Y21=Y(21);Y22=Y(22);Y23=Y(23);Y24=Y(24);Y25=Y(25);k1=k(1);k2=k(2);k3=k(3);k4=k(4);k5=k(5);k6=k(6);k7=k(7);k8=k(8);k9=k(9);k10=k(10);k11=k(11);k12=k(12);k13=k(13);k14=k(14);

dY1dt=-100*101.325/5/8.3145/773.15*(k1+k2+k3+k4+k5)*Y1/(Y1/644+Y2/230+Y3/115+Y4/52+Y5/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y1-Y2-Y3-Y4-Y5)/6)^(-0.385);

dY2dt=100*101.325/5/8.3145/773.15*(k1*Y1-(k8+k9+k10+k11)*Y2)/(Y1/644+Y2/230+Y3/115+Y4/52+Y5/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y1-Y2-Y3-Y4-Y5)/6)^(-0.385);

dY3dt=100*101.325/5/8.3145/773.15*(k2*Y1+k8*Y2-(k12+k13+k14)*Y3)/(Y1/644+Y2/230+Y3/115+Y4/52+Y5/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y1-Y2-Y3-Y4-Y5)/6)^(-0.385);

dY4dt=100*101.325/5/8.3145/773.15*(k3*Y1+k9*Y2+k12*Y3)/(Y1/644+Y2/230+Y3/115+Y4/52+Y5/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y1-Y2-Y3-Y4-Y5)/6)^(-0.385);

dY5dt=100*101.325/5/8.3145/773.15*(k4*Y1+k10*Y2+k13*Y3)/(Y1/644+Y2/230+Y3/115+Y4/52+Y5/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y1-Y2-Y3-Y4-Y5)/6)^(-0.385);

dY6dt=-100*101.325/10/8.3145/773.15*(k1+k2+k3+k4+k5)*Y6/(Y6/644+Y7/230+Y8/115+Y9/52+Y10/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y6-Y7-Y8-Y9-Y10)/6)^(-0.385);

dY7dt=100*101.325/10/8.3145/773.15*(k1*Y6-(k8+k9+k10+k11)*Y7)/(Y6/644+Y7/230+Y8/115+Y9/52+Y10/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y6-Y7-Y8-Y9-Y10)/6)^(-0.385);

dY8dt=100*101.325/10/8.3145/773.15*(k2*Y6+k8*Y7-(k12+k13+k14)*Y8)/(Y6/644+Y7/230+Y8/115+Y9/52+Y10/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y6-Y7-Y8-Y9-Y10)/6)^(-0.385);

dY9dt=100*101.325/10/8.3145/773.15*(k3*Y6+k9*Y7+k12*Y8)/(Y6/644+Y7/230+Y8/115+Y9/52+Y10/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y6-Y7-Y8-Y9-Y10)/6)^(-0.385);

dY10dt=100*101.325/10/8.3145/773.15*(k4*Y6+k10*Y7+k13*Y8)/(Y6/644+Y7/230+Y8/115+Y9/52+Y10/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y6-Y7-Y8-Y9-Y10)/6)^(-0.385);

dY11dt=-100*101.325/15/8.3145/773.15*(k1+k2+k3+k4+k5)*Y11/(Y11/644+Y12/230+Y13/115+Y14/52+Y15/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y11-Y12-Y13-Y14-Y15)/6)^(-0.385);

dY12dt=100*101.325/15/8.3145/773.15*(k1*Y11-(k8+k9+k10+k11)*Y12)/(Y11/644+Y12/230+Y13/115+Y14/52+Y15/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y11-Y12-Y13-Y14-Y15)/6)^(-0.385);

dY13dt=100*101.325/15/8.3145/773.15*(k2*Y11+k8*Y12-(k12+k13+k14)*Y13)/(Y11/644+Y12/230+Y13/115+Y14/52+Y15/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y11-Y12-Y13-Y14-Y15)/6)^(-0.385);

dY14dt=100*101.325/15/8.3145/773.15*(k3*Y11+k9*Y12+k12*Y13)/(Y11/644+Y12/230+Y13/115+Y14/52+Y15/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y11-Y12-Y13-Y14-Y15)/6)^(-0.385);

dY15dt=100*101.325/15/8.3145/773.15*(k4*Y11+k10*Y12+k13*Y13)/(Y11/644+Y12/230+Y13/115+Y14/52+Y15/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y11-Y12-Y13-Y14-Y15)/6)^(-0.385);

dY16dt=-100*101.325/20/8.3145/773.15*(k1+k2+k3+k4+k5)*Y16/(Y16/644+Y17/230+Y18/115+Y19/52+Y20/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y16-Y17-Y18-Y19-Y20)/6)^(-0.385);

dY17dt=100*101.325/20/8.3145/773.15*(k1*Y16-(k8+k9+k10+k11)*Y17)/(Y16/644+Y17/230+Y18/115+Y19/52+Y20/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y16-Y17-Y18-Y19-Y20)/6)^(-0.385);

dY18dt=100*101.325/20/8.3145/773.15*(k2*Y16+k8*Y17-(k12+k13+k14)*Y18)/(Y16/644+Y17/230+Y18/115+Y19/52+Y20/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y16-Y17-Y18-Y19-Y20)/6)^(-0.385);

dY19dt=100*101.325/20/8.3145/773.15*(k3*Y16+k9*Y17+k12*Y18)/(Y16/644+Y17/230+Y18/115+Y19/52+Y20/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y16-Y17-Y18-Y19-Y20)/6)^(-0.385);

dY20dt=100*101.325/20/8.3145/773.15*(k4*Y16+k10*Y17+k13*Y18)/(Y16/644+Y17/230+Y18/115+Y19/52+Y20/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y16-Y17-Y18-Y19-Y20)/6)^(-0.385);

dY21dt=-100*101.325/25/8.3145/773.15*(k1+k2+k3+k4+k5)*Y21/(Y21/644+Y22/230+Y23/115+Y24/52+Y25/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y21-Y22-Y23-Y24-Y25)/6)^(-0.385);

dY22dt=100*101.325/25/8.3145/773.15*(k1*Y21-(k8+k9+k10+k11)*Y22)/(Y21/644+Y22/230+Y23/115+Y24/52+Y25/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y21-Y22-Y23-Y24-Y25)/6)^(-0.385);

dY23dt=100*101.325/25/8.3145/773.15*(k2*Y21+k8*Y22-(k12+k13+k14)*Y23)/(Y21/644+Y22/230+Y23/115+Y24/52+Y25/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y21-Y22-Y23-Y24-Y25)/6)^(-0.385);

dY24dt=100*101.325/25/8.3145/773.15*(k3*Y21+k9*Y22+k12*Y23)/(Y21/644+Y22/230+Y23/115+Y24/52+Y25/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y21-Y22-Y23-Y24-Y25)/6)^(-0.385);

dY25dt=100*101.325/25/8.3145/773.15*(k4*Y21+k10*Y22+k13*Y23)/(Y21/644+Y22/230+Y23/115+Y24/52+Y25/16)/(1+0.45*k6)/(1+0.0989/6*k7)*(1+5.876*(100-Y21-Y22-Y23-Y24-Y25)/6)^(-0.385);

dYdt=[dY1dt;dY2dt;dY3dt;dY4dt;dY5dt;dY6dt;dY7dt;dY8dt;dY9dt;dY10dt;dY11dt;dY12dt;dY13dt;dY14dt;dY15dt;dY16dt;dY17dt;dY18dt;dY19dt;dY20dt;dY21dt;dY22dt;dY23dt;dY24dt;dY25dt];

运行结果如下：

Warning: Large-scale (trust region) method does not currently solve this type of problem,

switching to medium-scale (line search).

> In fmincon at 260

In k1_k14 at 16

Maximum number of function evaluations exceeded;

increase OPTIONS.MaxFunEvals.

使用函数fmincon()估计得到的参数值为:

k1 = 208.6494742904339100

k2 = 0.0000035022122300

k3 = 0.0093469868876100

k4 = 0.0176950218965225

k5 = 2919.9850027515131000

k6 =       280.852301

k7 = 33.2914927771165170

k8 = 0.0525562536448089

k9 = 0.1134399210673997

k10 = 0.0316754052171818

k11 = 0.6376342235568976

k12 = 946.3836064123888700

k13 = 16.5157190177285540

k14 = 146.0478939549901400

The sum of squares is: 544223.5729754469400000

Optimization terminated: relative function value

changing by less than OPTIONS.TolFun.

使用函数lsqnonlin()估计得到的参数值为:

k1 = 2511.6598737907384000

k2 = 1031.1081865545771000

k3 = 2192.9494255981772000

k4 = 2193.1454571674008000

k5 = 1967.9594464776890000

k6 =         0.115246

k7 = 4.0785589978752927

k8 = 9.2127248254471024

k9 = 13.8979528091439450

k10 = 8.3333906454971007

k11 = 0.4998132865030206

k12 = 2.0894143534582192

k13 = 0.5305261830822440

k14 = 0.4999974466911072

The sum of squares is: 863254.0758202144600000

方差达到6位数！！！！？？？？？是否与数据组数较少有关呢？也即产生了多解导致的方差太大？

Xsim=Xsim'是自己加上去的，不然程序报错，矩阵大小不一致！