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];             %參數(shù)初值
lb=[0 0 0 0 0 0 0 0];  ub = [1 1 1 1 1 1 1 1]*1e1;  %lb、ub：參數(shù)下限和上限
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);            

%使用函數(shù)fmincon()進(jìn)行參數(shù)估計(jì)
[k,fval,flag]=fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],Y0,yexp);
fprintf('\n 使用函數(shù)fmincon()估計(jì)得到的參數(shù)值為:\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的計(jì)算值','Y1的實(shí)驗(yàn)值','Y2的計(jì)算值','Y2的實(shí)驗(yàn)值','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];