function ode_fit
format long
clear all
clc

x0 = [6.7 4.86 3.17];
k0 = [1.5 1.5 1.5 6.8 2/3 0.1 0.1 0.01];   
lb = -[1 1 1 1 1 1 1 1 ]*1e9;
ub = [1 1 1 1 1 1 1 1 ]*1e9;

data=...
    [

88        6.7        3.17        4.86
92        7.5        2.87        5.17                                   
96        8.4        2.33        5.72                 
100        9.1        1.87        6.18
104        9.3        1.36        6.73               
108        9.1        0.95        7.13
112        8.6        0.52        7.56
];
tspan = [data(:,1)'-88];
yexp = [data(2:end,2) data(2:end,4) data(2:end,3)];

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.9f \n',k(1))
fprintf('\tk2 = %.9f \n',k(2))
fprintf('\tk3 = %.9f \n',k(3))
fprintf('\tk4 = %.9f \n',k(4))
fprintf('\tk5 = %.9f \n',k(5))
fprintf('\tk6 = %.9f \n',k(6))
fprintf('\tk7 = %.9f \n',k(7))
fprintf('\tk8 = %.9f \n',k(8))

figure(1)

ts=0:((max(tspan)-min(tspan))/100):max(tspan);
[ts ys] = ode45(@KineticsEqs,ts,x0,[],k);
yy = [data(:,2) data(:,4) data(:,3)];
figure(1)
plot(ts,ys(:,1),'b',tspan,yy(:,1),'bo',ts,ys(:,2),'r',tspan,yy(:,2),'ro',ts,ys(:,3),'k',tspan,yy(:,3),'ko'),
legend('X的计算值','X的实验值','P的计算值','P的实验值','S的计算值','S的实验值','Location','best');



function f = ObjFunc(k,tspan,x0,yexp)           % 目标函数
[t Xsim] = ode45(@KineticsEqs,tspan,x0,[],k);
Xsim1=Xsim(:,1);
Xsim2=Xsim(:,2);
Xsim3=Xsim(:,3);
ysim(:,1) = Xsim1(2:end);
ysim(:,2) = Xsim2(2:end);
ysim(:,3) = Xsim3(2:end);

size(ysim(:,1));
size(ysim(:,2));
size(yexp(:,1));
size(yexp(:,2));

f = [(ysim(:,1)-yexp(:,1)) (ysim(:,2)-yexp(:,2)) (ysim(:,3)-yexp(:,3))];

function dCdt = KineticsEqs(t,C,k)              % ODE模型方程
umax=k(1);Ks=k(2);Ki=k(3);Pmax=k(4);a=k(5);Yxs=k(6);Yps=k(7);m=k(8);
X=C(1);P=C(2);S=C(3);

dXdt=umax*(S./(S+Ks)).*(1./(1+S/Ki)).*(1-P/Pmax).*X;
dPdt=a*dXdt;
dSdt=-Yxs*dXdt-Yps*dPdt-m*X;

dCdt = [dXdt; dPdt;dSdt];