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

4Â¥: Originally posted by ygqsgc at 2021-07-18 16:37:00
лл£¬µ«ÊÇÎÒÐÞ¸ÄÕâÀïºÃÏñûʲô±ä»¯¡£ÎÒ½ñÌìÓÖ·¢ÏÖ³õÖµ±ä»¯»áºÜ´óµÄÓ°Ïì½á¹û£¬ÇëÎÊÄúÓÐʲô½¨ÒéÂð£¿...

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)
¼ì²éÒ»ÏÂÕâÀÓÈÆäÊǵÚÒ»ÐÐ