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)
检查一下这里，尤其是第一行