ode15s error（Genetic algorithm for differential equation parameter identification）

5 次查看（过去 30 天）

Bowen Yu 2021-1-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/716638-ode15s-error-genetic-algorithm-for-differential-equation-parameter-identification

评论： Star Strider 2021-6-20

Hello. I Truly appreciate if somebody check my code. I want to use the experimental data to identify the unknown parameters of the differential equation model.

Here’s my code：

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

Fitness function：

function [Sfit] = my_fitnessfun(a)
tspan=[0 10 20 30 40 50 60 70 80 90 100 110 120 130];
% Value of experiment
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];
% Initial value
y0(1)=Sreal(1);
y0(2)=Ereal(1);
y0(3)=Xreal(1);
y0(4)=Greal(1);
% To solve the differential equation, we introduce the differential equation function@myfun  
options = odeset('RelTol',1e-3,'AbsTol',1e-6);
[t,ycal]=ode15s(@(t,y) myfun(t,y,a),tspan,y0,options);
% 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

Differential equation：

function dydt = myfun(t, y, a)
%18 parameters
Kh = a(1);%rate constant
Km = a(2);%Michaelis constant
KG = a(3);%glucose inhibition constant
Kstarch = a(4);%starch inhibition constant
qm = a(5);%maximum specificproductionrate
Ksp = a(6);%saturation production constant
Kssp = a(7);%substrate inhibition constant
Kex = a(8);%product inhibition constant
mium= a(9);%maximum specific growth rate of biomass
Ks = a(10);%saturation growth constant
Kss = a(11);%substrate inhibition constant
Kxx = a(12);%product inhibition constant
Kd = a(13);%cell death constant
YXG = a(14);%yield coefficient of biomass growth
YEG = a(15);% yield coefficient of ethanol
% beta = a(16);%enzyme degradation rate
% Kenz = a(17);%enzyme inhibition constant 
% Enzm = a(18);%maximum enzyme concentration
%S, E, X, G, Enz were starch, ethanol, cell, enzyme concentration
S=y(1);
E=y(2);
X=y(3);
G=y(4);
% Enz=y(5);
Enz=1;%Without the enzyme concentration
%Starch consumption rate
Rs = Kh*Enz/(Km*(1 + G/KG) + (S^2)/Kstarch + S);
%Cell specific growth rate
miu = mium*G*exp(-E/Kxx)/(Ks + G + (G^2)/Kss);
%Ethanol production rate
q = qm*G*exp(-E/Kex)/(Ksp + G + (G^2)/Kssp);
% %Rate of enzyme synthesis
% Renz = (mium + beta)*Enzm*S/(Kenz + S);
dydt = zeros(4,1);
%Starch utilization
dydt(1) = -Rs*S;
%Ethanol kinetics
dydt(2) = miu*X - Kd*X;
%Growth kinetics
dydt(3) = q*X;
%Glucose utilization
dydt(4) = 1.111*Rs*S - (1/YXG)*dydt(2) - (1/YEG)*dydt(3);
% %Enzyme kinetics
% dydt(5) = Renz - (miu + beta)*Enz;
end

Error after running code：

Index in position 1 exceeds array bounds (must not exceed 1).

Error in my_fitnessfun (line 25)

ff(l)=((ycal(l,1)-Sreal(l))/max(Sreal))^2 + ...

Error in createAnonymousFcn>@(x)fcn(x,FcnArgs{:}) (line 11)

fcn_handle = @(x) fcn(x,FcnArgs{:});

Error in makeState (line 52)

firstMemberScore = FitnessFcn(state.Population(initScoreProvided+1,:));

Error in gaunc (line 40)

state = makeState(GenomeLength,FitnessFcn,Iterate,output.problemtype,options);

Error in ga (line 398)

[x,fval,exitFlag,output,population,scores] = gaunc(FitnessFcn,nvars, ...

Caused by:

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

The position is [ error my_fitnessfun (line 25)] , so running at the break point on this line, the size of ycal before the error is [1,1] , and the later for loop needs to refer to ycal (L, 1) , the size of L ranges from 1 to 14, it would be out of the range of YCAL’s index (the numbers here are just an example, and the numbers will be slightly different from one error to another) , ycal is calculated by ode15s, but before that Ode15s warned that they could not compute the integral, so check the YCAL and T output of Ode15s and find that the T dimension of the output is different from tspan due to numerical calculation problems. Fewer t will be returned when the ODE fails halfway through.

I found out what was wrong with my code, but I didn’t know how to fix it. I Truly appreciate if somebody check my code.

3 个评论
显示 1更早的评论隐藏 1更早的评论

Alex Sha 2021-1-23

Hi, Yu, refer to the results below please:

Weighted Root of Mean Square Error (RMSE): 0.403256606110099
Weighted Sum of Squared Residual: 8.45602629931463
Correlation Coef. (R): 0.995510242685548
R-Square: 0.991040643291839
Parameter          	Best Estimate
--------------------	-------------
kh	-0.420591358564488
km	-5.93463278095738
kg	0.0472576188654033
kstarch	-21.5825845231028
mium	-36.0337479065289
kxx	3.41142942705453
ks	-19.039415556528
kss	-0.00591853503606009
kd	0.00791615631803505
qm	-0.00726897451832187
kex	5.00111221131158
ksp	-0.123904200107017
kssp	-0.472693054329812
yxg	0.156258343283814
yeg	1.04741098630007

Bowen Yu 2021-1-25

hello！But how did you get that result

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Star Strider 2021-1-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/716638-ode15s-error-genetic-algorithm-for-differential-equation-parameter-identification#answer_597928

I cannot follow your code.

See Ode system solution with unknown constant for one example of an approach that generally works.

25 个评论
显示 23更早的评论隐藏 23更早的评论

Star Strider 2021-2-18

This line:

ftns = @(theta) norm(y-solver(theta,t));

creates the fitness function. It takes the norm of the result of subtracting the output of the objective function (containing the differential equation and the solution), and returns the scalar result for the specified values of the coefficient vector ‘theta’.

There are several possible reasons of the problems you are encountering. Most likely is that the parameters are not close to the correct magnitude, so the differential equation encounters a singularity (±Inf) result, and stops integrating.

I ran ‘bowen_15.m’ for this analysis.

There could be something wrong with the system itself, since the data do not appear to even closely approximate the results of the integration.

I was able to prevent the problems with the integration not completing (at least most of the time) by changing the population size and scaling the initial poulation matrix to a lower value:

PopSz = 25;
Parms = 15;
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-6, 'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);

I also generally let the optimisation routine choose the initial values , so within ‘solver’, I set:

y0 = a(16:19);

instead of using the pre-determined initial values, and in the calling code set:

Parms = 19;

to accommodate that change.

Experiment with that to see if it improves things.

I also added:

fprintf('\nBest Fitness = %11.6f\n',fval)

just before printing the parameters. This helps to keep track of it.

So I encourage you to search for problems in the differential equation system you coded, since it apparently is not actually able to model the data you have. Over several runs, most of the results should approximate the data. That is not the situation here.

Bowen Yu 2021-6-17

I’m sorry. My problem isn’t solved. I’m still gonna need your help.

I referenced an article and tested my code against the data in the article.

Same as before. When fitting a single differential equation, it works very well, but when fitting a differential equation, it doesn’t work, or sometimes it does, but it doesn’t work.(Code file uploaded)

This is the single differential equation mentioned above：

This is the differential equation in the code：

t is actually time, S, G, X, E is the experimental data. The rest are constants to be identified.

%Experimental data
t = [0 1 2 5 7 9 13 20 24 30 36 38 40 44 50 60 64]';
S = [100 95 90 85 78 70 55 40 30 20 11 10.5 10 9.5 7 2 0]';
E = [0 0.1 1 1.5 2 2.5 7 12.5 15 20 22.5 24.5 25 25.2 29 30 30.1]';
X = [0.609 0.65 0.655 0.7 0.75 0.8 1.15 1.5 1.8 2.2 2.3 2.32 2.35 2.32 2.4 2.5 2.4]';
G = [0 2 6 10 14 16 15 10 7 4 3.5 3 2.5 2 1.5 0.5 0.2]';

This is what I found when I tested it with a single differential equation：（It’s very good, but when you have multiple equations, it doesn’t work）

I don’t know what’s wrong with my code, but I could really need your help.

Thank you so much for helping me.

Alex Sha 2021-6-20

The code for obtaining above results in 1stOpt are as below:

WeightedReg = 3;
Variable t,S,E,X,G;
ODEFunction
S'=-Vm*S/(Km+S);
G'=1.111*Vm*S/(Km+S)-(1/Yps)*E'-(1/Yxs)*X'-Ms*X;
X'=(Mium*G/(Ks+G))*X;
E'=(Qm*G/(Ks+G))*X;
Data;
t = [0 1 2 5 7 9 13 20 24 30 36 38 40 44 50 60 64];
S = [100 95 90 85 78 70 55 40 30 20 11 10.5 10 9.5 7 2 0];
E = [0 0.1 1 1.5 2 2.5 7 12.5 15 20 22.5 24.5 25 25.2 29 30 30.1];
X = [0.609 0.65 0.655 0.7 0.75 0.8 1.15 1.5 1.8 2.2 2.3 2.32 2.35 2.32 2.4 2.5 2.4];
G = [0 2 6 10 14 16 15 10 7 4 3.5 3 2.5 2 1.5 0.5 0.2];