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

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

4楼: Originally posted by ygqsgc at 2021-07-18 16:37:00
谢谢，但是我修改这里好像没什么变化。我今天又发现初值变化会很大的影响结果，请问您有什么建议吗？...

5楼: Originally posted by 独孤神宇 at 2021-07-18 17:32:29
这个需要多次尝试了，可以将结果作为初始值代入反复计算。

可以参考三楼的结果，将其作为 初值代入计算...

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