function Kinetics2
clear all; clc
format long
k0=[0 0 0 0 0 0 0 0];
tspan=[0 2 4 6 8 10 15 20 25 30 40 50 60 75 90 105 120];             %²ÎÊý³õÖµ
lb=[0 0 0 0 0 0 0 0];  ub = [1 1 1 1 1 1 1 1]*1e1;  %lb¡¢ub£º²ÎÊýÏÂÏÞºÍÉÏÏÞ
Y0=[0 0];
Kinetics1=[0        0        0
2        0.101406296        0.025750304
4        0.148139381        0.077208075
6        0.171971874        0.097848253
8        0.200101464        0.131543106
10        0.223088343        0.149258466
15        0.268349285        0.206242074
20        0.319194899        0.249534539
25        0.352139308        0.280899049
30        0.373849138        0.309100951
40        0.4145939        0.348550324
50        0.442947836        0.381294016
60        0.460584916        0.395704095
75        0.483226647        0.421979897
90        0.499672966        0.430587898
105        0.50291736        0.452012783
120        0.518449036        0.452107899];
yexp=Kinetics1(1:17,2:3);            

%ʹÓú¯Êýfmincon()½øÐвÎÊý¹À¼Æ
[k,fval,flag]=fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],Y0,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('\tk6=%.4f\n',k(6)),fprintf('\tk7=%.4f\n',k(7)),fprintf('\tk8=%.4f\n',k(8));
fprintf(' The sum of the squares is:%.1e\n\n', fval);
k_fmincon=k;

ts=0:1:max(tspan);
[ts ys] = ode45(@func,ts,Y0,[],k);
yy = [Kinetics1(:,2) Kinetics1(:,3)];
plot(ts,ys(:,1),'b',tspan,yy(:,1),'bo',ts,ys(:,2),'r',tspan,yy(:,2),'ro'),
legend('Y1µÄ¼ÆËãÖµ','Y1µÄʵÑéÖµ','Y2µÄ¼ÆËãÖµ','Y2µÄʵÑéÖµ','Location','best');

function f=ObjFunc4Fmincon(k,Y0,yexp)
tspan=[0 2 4 6 8 10 15 20 25 30 40 50 60 75 90 105 120];
[t,Y]=ode45(@func,tspan,Y0,[],k);
y(:,1)=Y(:,1);y(:,2:3)=Y(:,1:2);
f=sum((y(:,2)-yexp(:,1)).^2)+sum((y(:,3)-yexp(:,2)).^2)

function dYdt=func(t,Y,k)
k1=k(1);k2=k(2);k3=k(3);k4=k(4);k5=k(5);k6=k(6);k7=k(7);k8=k(8);
f1=(k1.*k2.*k3.*k4.*(1-Y(1)).*6.25.*(1-Y(2)).*3.125-k5.*k6.*k7.*k8.*Y(1).*6.25.*Y(2).*3.125)./(((k1.*k2.*k7+k1.*k2.*k4).*(1-Y(1)).*6.25+(k5.*k3.*k4+k2.*k3.*k4).*(1-Y(2)).*3.125+(k5.*k6.*k7+k1.*k6.*k7).*Y(1).*6.25+(k2.*k6.*k8+k5.*k7.*k8).*Y(2).*3.125+(k1.*k3.*k4+k1.*k2.*k3).*(1-Y(1)).*6.25.*(1-Y(2)).*3.125+(k1.*k6.*k7+k5.*k6.*k4).*(1-Y(1)).*6.25.*Y(1).*6.25+(k2.*k4.*k8+k5.*k3.*k8).*(1-Y(2)).*3.125+(k6.*k7.*k8+k5.*k6.*k8).*Y(1).*6.25.*Y(2).*3.125).*6.125);
f2=(k1.*k2.*k3.*k4.*(1-Y(1)).*6.25.*(1-Y(2)).*3.125-k5.*k6.*k7.*k8.*Y(1).*6.25.*Y(2).*3.125)./(((k1.*k2.*k7+k1.*k2.*k4).*(1-Y(1)).*6.25+(k5.*k3.*k4+k2.*k3.*k4).*(1-Y(2)).*3.125+(k5.*k6.*k7+k1.*k6.*k7).*Y(1).*6.25+(k2.*k6.*k8+k5.*k7.*k8).*Y(2).*3.125+(k1.*k3.*k4+k1.*k2.*k3).*(1-Y(1)).*6.25.*(1-Y(2)).*3.125+(k1.*k6.*k7+k5.*k6.*k4).*(1-Y(1)).*6.25.*Y(1).*6.25+(k2.*k4.*k8+k5.*k3.*k8).*(1-Y(2)).*3.125+(k6.*k7.*k8+k5.*k6.*k8).*Y(1).*6.25.*Y(2).*3.125).*3.125);
dYdt=[f1;f2];

3Â¥: Originally posted by dingd at 2021-07-18 10:50:32
ÓÃ1stOptÊÔÁËÏ£¬ËƺõÊǶà½â°É£º

1£º
¾ù·½²î(RMSE): 0.133593203954135
²Ð²îƽ·½ºÍ(SSR): 0.571108612567395
Ïà¹ØϵÊý(R): 0.999611248630951
Ïà¹ØϵÊý֮ƽ·½(R^2): 0.999222648389528
ÐÞÕýRƽ·½(Adj. R^2) ...

2Â¥: Originally posted by ¶À¹ÂÉñÓî at 2021-07-17 23:47:54
y(:,1)=Y(:,1);y(:,2:3)=Y(:,1:2);
f=sum((y(:,2)-yexp(:,1)).^2)+sum((y(:,3)-yexp(:,2)).^2)
¼ì²éÒ»ÏÂÕâÀÓÈÆäÊǵÚÒ»ÐÐ