Genetic Algorithm for differential equation parameter identification

To estimate the differential equation parameters a 1, a 2, a 3. Here's what I wrote, but didn't get the expected result (I differential equation the initial a value 1,2,4 and got the s value as the actual experimental data, so the estimated result should still be 1,2,4) , i hope the boss can help me find the wrong part

function myproject
clear;
clc;

%ga
options = optimoptions('ga','MaxGenerations',1000,'MaxStallGenerations',300,...
     'MaxStallTime',50,'FunctionTolerance',1e-12,'ConstraintTolerance',1e-12);
[a0,fval,exitflag,output]=ga(@my_fitnessfun3,3,options);
fprintf('\n\nEstimation value of genetic algorithm:\n');
disp(a0);
end

function [Sfit] = my_fitnessfun3(a)%Fitness function

% tspan
tspan=[0 14 28 42 56 70 84 98 112 126 140]';

% S True value of experiment
Sreal=[30 31.5772657565012 33.0963189238513 34.5639532181971 35.9811192405661 37.3514763949577 38.6783813122688 39.9655204641298 41.2159922182303 42.4329666567832 43.6190705086730]';

% S Initial value
S0=30;

% Give the parameter an initial value
a=[1 2 4];

% To solve the differential equation, we introduce the differential equation function@myfun3 and assign a value to it
[t,Scal]=ode15s(@myfun2,tspan,S0,[],a);

% minf(l)
n=length(Sreal);

for l=1:n
   ff(l)=(Scal(l)-Sreal(l))^2;
end
Sfit=sum(ff(l));
end

function [dSdt] = myfun2(t,S,a)%DIFFERENTIAL EQUATION

% Three parametersa(1),a(2),a(3)
mjumax = a(1);
Ks = a(2);
K1 = a(3);

dSdt = (mjumax*S)/(S+Ks+S^2/K1);
end

Answer 1

[Alan Weiss]

Firstly, the code did not work as you gave it. The function you labeled as myfun2 should have been labeled myfun3.

Secondly, your fitness function code my_fitnessfun3 overwrites the input variable a with the value [1 2 4], so every input value a is ignored. Just comment out the line setting a.

Thirdly, ODE solvers do not behave well when given bad inputs. Set lower bounds of 0 and upper bounds of, say, 20 on all variables. Then things run sensibly.

options = optimoptions('ga','MaxGenerations',1000,'MaxStallGenerations',300,...
     'MaxStallTime',50,'FunctionTolerance',1e-12,'ConstraintTolerance',1e-12);
[a0,fval,exitflag,output]=ga(@my_fitnessfun3,3,[],[],[],[],[0 0 0],[20 20 20],[],options)

Finally, you get much a more accurate answer much more quickly by using the correct solver for the problem, lsqnonlin or lsqcurvefit, rather than ga. See Fit an Ordinary Differential Equation.

Alan Weiss

MATLAB mathematical toolbox documentation

[Bowen Yu]

Thank you for your answer. But why do I get different results every time I run, should I not use GA to solve this kind of problem?

[Alan Weiss]

As I said, "you get much a more accurate answer much more quickly by using the correct solver for the problem, lsqnonlin or lsqcurvefit, rather than ga. See Fit an Ordinary Differential Equation." If you change to lsqnonlin or lsqcurvefit, you need to modify your objective function to not return the sum of squares of the Scal - Sreal vector, but to return the entire vector.

Alan Weiss

MATLAB mathematical toolbox documentation

[Alan Weiss]

One more thing. I believe that a large number of parameters a solve the problem. This means that there is no unique solution; many parameterizations lead to the same ODE solution.

Alan Weiss

MATLAB mathematical toolbox documentation

[Bowen Yu]

Thank you so much for your help！