GA辨识微分方程参数问题

clear;
clc;
%ga
options = optimoptions('ga','MaxGenerations',200,'MaxStallGenerations',50,...
     'MaxStallTime',inf,'FunctionTolerance',1e-12,'ConstraintTolerance',1e-6);
[a0,fval,exitflag,output]=ga(@my_fitnessfun,15,[],[],[],[],...
    [],[],[],options)


运行后报错：

警告: 在 t=2.888835e+00 处失败。在时间 t 处，步长大小必须降至所允许的最小值(7.105427e-15)以下，才能达到积分容差要求。
位置 1 处的索引超出数组边界(不能超出 1)。

出错 my_fitnessfun (line 25)
    ff(l)=((ycal(l,1)-Sreal(l))/max(Sreal))^2 + ...

出错 createAnonymousFcn>@(x)fcn(x,FcnArgs{:}) (line 11)
fcn_handle = @(x) fcn(x,FcnArgs{:});

出错 makeState (line 52)
            firstMemberScore = FitnessFcn(state.Population(initScoreProvided+1,:));

出错 gaunc (line 40)
state = makeState(GenomeLength,FitnessFcn,Iterate,output.problemtype,options);

出错 ga (line 398)
            [x,fval,exitFlag,output,population,scores] = gaunc(FitnessFcn,nvars, ...

原因:
    Failure in initial user-supplied fitness function evaluation. GA cannot continue.





以下是适应度函数
function [Sfit] = my_fitnessfun(a)

% t取值区间
tspan=[0 10 20 30 40 50 60 70 80 90 100 110 120 130];

% 实验真实值
Sreal=[28.5 27.1 23.2 20.0 16.6 14.3 11.5 9.2 7.9 6.9 5.2 3.6 2.0 1.5];
Ereal=[0 0.24 1.06 2.20 3.10 4.06 5.26 6.31 7.23 7.82 7.94 8.09 8.18 8.20];
Xreal=[0.14 0.24 0.54 1.05 1.34 1.64 1.97 2.31 2.82 3.04 3.13 3.22 3.12 3.09];
Greal=[0 0.423 0.184 0.174 0.153 0.146 0.136 0.133 0.127 0.123 0.117 0.108 0.103 0.095];

% 初始值
S0(1)=Sreal(1);
S0(2)=Ereal(1);
S0(3)=Xreal(1);
S0(4)=Greal(1);

% 解微分方程，引入微分方程函数@myfun
options = odeset('RelTol',1e-3,'AbsTol',1e-6);
[t,ycal]=ode15s(@myfunn,tspan,S0,options,a);

% minf(l)
n=length(tspan);
for l=1:n
    ff(l)=((ycal(l,1)-Sreal(l))/max(Sreal))^2 + ...
        ((ycal(l,2)-Ereal(l))/max(Ereal))^2 + ...
        ((ycal(l,3)-Xreal(l))/max(Xreal))^2 + ...
        ((ycal(l,4)-Greal(l))/max(Greal))^2;
end
Sfit=sum(ff(l));
end
以下是微分方程
function [dydt] = myfunn(t,y,a)

%18个参数
Kh = a(1);%速率常数
Km = a(2);%米氏常数
KG = a(3);%葡萄糖抑制常数
Kstarch = a(4);%淀粉抑制常数
qm = a(5);%最大产率
Ksp = a(6);%饱和生产常数
Kssp = a(7);%底物抑制常数
Kex = a(8);%产物抑制常数
mium= a(9);%最大比生长速率
Ks = a(10);%饱和生长常数
Kss = a(11);%底物抑制常数
Kxx = a(12);%产物抑制常数
Kd = a(13);%细胞死亡常数
YXG = a(14);%生物量生长产率系数
YEG = a(15);%乙醇产率系数
% beta = a(16);%酶降解率
% Kenz = a(17);%酶抑制常数
% Enzm = a(18);%最大酶浓度

%S,E,X,G,Enz分别为淀粉，乙醇，细胞，酶的浓度
S=y(1);
E=y(2);
X=y(3);
G=y(4);
% Enz=y(5);
Enz=1;

dydt = zeros(4,1);
%淀粉利用
dydt(1) = -Kh*Enz*S/(Km*(1 + G/KG) + S^2/Kstarch + S);

%乙醇动力学
dydt(2) = qm*G*X*exp(-E/Kex)/(Ksp + G + G^2/Kssp);

%细胞生长动力学
dydt(3) = mium*G*exp(-E/Kxx)*X/(Ks + G + G^2/Kss) - Kd*X;

%葡萄糖利用
dydt(4) = 1.111*Kh*Enz/(Km*(1 + G/KG) + S^2/Kstarch + S) - ...
    (1/YXG)*mium*G*exp(-E/Kxx)/(Ks + G + G^2/Kss) - ...
    (1/YEG)*qm*G*X*exp(-E/Kex)/(Ksp + G + G^2/Kssp);

%酶动力学
% dydt(5) = (mium + beta)*Enzm*S/(Kenz + S) - ...
%     (mium*G*exp(-E/Kxx)/(Ks + G + G^2/Kss) + beta)*Enz;
end



真的是弄的我焦头烂额。。。有没有大佬教教我小菜鸟，有偿也可以。。。拜托了