ode15s error（Genetic algorithm for differential equation parameter identification）

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.

[Bowen Yu]

This is my code file

[Alex Sha]

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]

hello！But how did you get that result

Answer 1

[Star Strider]

I cannot follow your code.

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

[Bowen Yu]

Hello ,thank you for your help! Could you please post your ga code for me .^-^

[Star Strider]

My pleasure!

I just now attached it to my last Comment to that Answer, since it belongs there.

[Bowen Yu]

Sorry to reply you so late because of something.

Thank you so much for your help!

I imitated your code and found that when there are few constants (bowen_2.m), the effect is very good. But when the differential equation becomes more complicated and the constant becomes more (bowen_15.m), the effect is very bad (use the resultant constant to graph the comparison with the experimental data). Moreover, the code does not run successfully every time, and sometimes it will report an error. The content of the error is the same as the code I wrote before.

You said before that you cannot follow my code, so I revised my code again. （aaa.m）

I imitated your code is also attached in the comments.（bowen_2.m, bowen_15.m）

I am very grateful for your help and hope you can help me again. >_<

[Bowen Yu]

Sorry. There is a part of your code I don’t know what it means. But I still imitated

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

[Bowen Yu]

The content of aaa.m and bowen_15.m errors are almost the same, but aaa.m will report an error every time, bowen_15.m is not every time, sometimes it can run successfully, but the result is not good.

[Star Strider]

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]

Thank you again for your help, the problem may be in the differential equation model as you said, it may be that the model and the data do not match well

[Star Strider]

As always, my pleasure!

I agree. I will continue to help once you get the problem with the differential equation system sorted. I cannot help with that, since I do not know what model you are using or the system you are estimating.

[Bowen Yu]

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.

[Star Strider]

There are errors in ‘dydt(2)’ that I corrected. I could not have caught those without ‘Table 3’, so I appreciate your including it.

This appears to be correct:

[~,yv] = ode45(@Difeq,t,y0);
    function dydt = Difeq(~,y)
    %Differential equation
    dydt = zeros(4,1);

    dydt(1) = -Vm * y(1) / (Km + y(1));

    dydt(2) = 1.111 * Vm * y(1) / (Km + y(1)) - (1/YXS) * dydt(4) - (1/YPS) * dydt(3) - ms * y(3);

    dydt(3) = (mium * y(2)  / (Ks + y(2))) * y(3);

    dydt(4) = (qm * y(2) / (Ks + y(2))) * y(3); 
    end
Y = yv;
end

In ‘solver’, I also defined:

y0 = [100 0 0.6 0];

However while I can get it to fit ‘S’ in almost every attempt, I cannot get it to fit anything else even closely. The model may not be appropriate for the experimental results, the units are not compatible (the data need to be scaled appropriately), or there could be other problems. Any other information would be helpful.

[Bowen Yu]

I really appreciate your help.

But I don’t think it’s a problem with the models or the data. The data and the model are all from the same article, but it is fitted using other methods. There must be something wrong with the code.

Thanks again for your kind help！

[Star Strider]

As always, my pleasure!

The posted model did appear to be different from the expressions in the table, so I corrected the model.

What parameter values does the article provide? What are the units of the concentrations?

Can you post a PDF of the article?

.

[Bowen Yu]

This is the article. Anyway, thank you very much!

[Bowen Yu]

There’s another article. I know from this article that genetic algorithms can be used to solve this problem, but he did not give a detailed

[Alex Sha]

Yu, how about the results below:

Root of Mean Square Error (RMSE): 1.11045737774729
Sum of Squared Residual: 78.9193976187769
Correlation Coef. (R): 0.992344112924485
R-Square: 0.984746838455883

Parameter          	Best Estimate
--------------------	-------------
vm	5.76284101843475
km	59.4316899561422
yps	0.206733436735961
yxs	0.209532706153286
ms	-0.519655541510672
mium	0.0744294531262563
ks	13.9141483401913
qm	1.64155965922936

[Bowen Yu]

Wow!!! Thank you so much! I finally figured out what the problem is. I got the same result, but why is X not so good.

[Alex Sha]

Hi, it is because the magnitude of the X data values are much less than others. It depends on the objective function selected, for example, If considering the weight factor，the results will be as below. The comparative pictures looks much better than before, especially for X, however, note the value of objective function "Sum of Squared Residual", the previous (78.9193976187769) is less than currently one （174.349592845946）

Weighted Root of Mean Square Error (RMSE): 1.6505188239514
Weighted Sum of Squared Residual: 174.349592845946
Correlation Coef. (R): 0.995148063403394
R-Square: 0.990319668095526

Parameter          	Best Estimate
--------------------	-------------
vm	3.96863007888113
km	27.6612817553871
yps	0.385331357873195
yxs	0.0450150770244026
ms	-0.0681061360649869
mium	0.21595543006129
ks	44.3571324582003
qm	3.19692235754256

[Star Strider]

I tried several different MATLAB solvers over several hours, the most recent being surrogateopt, and still cannot get the system of differential equations to converge on the correct result.

[Alex Sha]

Hi，@Star Strider, based on my own experience and many comparative tests，for global optimization, at least so far, no any other optimization solver can compare with 1stOpt. I would like to give you some test questions if you want, those questions are much hard than some classical problems, for example, the six-hump camel back function or Rosenbrock's function.

[Bowen Yu]

Hello,@Star Strider,My previous mistake was that the second equation should be written as follows. Without changing the other code, you can get almost the same result as @Alex Sha did after just doing this.

Thanks again to @Star Strider and @Alex Sha for your enthusiastic help, Finally, the problem of my master's thesis was finally solved. My previous major was bioengineering, so I am a beginner in MATLAB. Thanks for your help ！^_^

[Star Strider]

As always, my pleasure!

Thank you!

With those changes, it converged quickly and gave appropriate results! Since ga always works in situations like this, I was wondering why it failed here. With the correct system, it succeeded!

Where was that expression in the papers? I did not see it anywhere.

The correct system:

    function dydt = Difeq(~,y)
    %Differential equation
    % S = y(1),  G = y(2),  X = y(3),  E = y(4)
    dydt = zeros(4,1);

    dydt(1) = -Vm * y(1) / (Km + y(1));

    dydt(2) = 1.111 * Vm * y(1) / (Km + y(1)) - (1/YXS) * (mium * y(2))/(Ks+y(2)) * y(3) - (1/YPS) * (qm * y(2))/(Ks+y(2)) * y(3) - ms * y(3);
    

    dydt(3) = mium * y(2) * y(3) / (Ks + y(2));

    dydt(4) = qm * y(2) * y(3) / (Ks + y(2)); 
    end

.

[Bowen Yu]

There is no such expression in the article. I only know when I see Alex Sha's answer.

[Star Strider]

Thank you. He must either have derived it himself, or gotten it from another source.

[Alex Sha]

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

Press "View ODE Function", the follow equation will be generated automatically:

[Star Strider]

Fascinating!