How do I implement a genetic algorithm for parameter fitting to an already existing ODE system>
显示 更早的评论
So I am working on expanding an already working ODE system. I have my rate equations and fluxes for three additional differential equations. I would like to estimate 8 unique constants in the 3 rate equations which is all building upon a 20 equation ODE system. I also have vectorized how I would like my new 3 equations to behave for a specific input. How can I implement the genetic algorithm to focus on parameter fitting for this. I can't reveal specifically each portion of my code but i can put some filler code to demonstrate my issue. thanks in advance!
%raw data of how i want the model to behave
minutes = [10, 15, 20, 40, 60, 80]';
neweq1= [60,20,10,6,5,4.9 ]';
neweq2=[40,90,95,99,99.5]';
neweq3= [60,70,20,19.5,19]';
%fitfun= norm(fun(t)-expected value)
coeff= ga(fitfun,8);
function o= fitness (x,t,y)
y0=[200 0 .17 209.9 94.83 483.73 1.51 14.76 0 200 184.7 15.31 750 0 0 0 300 0 3000 0 100 0 0];
tspan=[0 100];
[t,y]=ode45 (@(t,y) fun, tspan, y0)
function dy=fun(t,y,input)
dy=zeros(23,0);
%20 rate equations
%my new 3 equations
%x is an array populated with the parameters i want to estimate and is a part of my new 3 equations
end
end
采纳的回答
Here is some prototype code I developed a while ago that you can adapt to your system —
t=[0.1
0.2
0.4
0.6
0.8
1
1.5
2
3
4
5
6];
c=[0.902 0.06997 0.02463 0.00218
0.8072 0.1353 0.0482 0.008192
0.6757 0.2123 0.0864 0.0289
0.5569 0.2789 0.1063 0.06233
0.4297 0.3292 0.1476 0.09756
0.3774 0.3457 0.1485 0.1255
0.2149 0.3486 0.1821 0.2526
0.141 0.3254 0.194 0.3401
0.04921 0.2445 0.1742 0.5277
0.0178 0.1728 0.1732 0.6323
0.006431 0.1091 0.1137 0.7702
0.002595 0.08301 0.08224 0.835];
ftns = @(theta) norm(c-kinetics(theta,t));
PopSz = 50;
Parms = 10;
optsAns = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10); % Options Structure For 'Answers' Problems
% opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
% optshf = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'HybridFcn',@fmincon, 'PlotFcn',@gaplotbestf, 'PlotInterval',1); % With 'HybridFcn'
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],optsAns);
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)
fprintf('Fitness value at convergence = %.4f\nGenerations \t\t = %d\n\n',fval,output.generations)
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
CV(k1,:) = hd(k1).Color;
hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','N')
function C=kinetics(theta,t)
c0 = theta(7:10);
% c0=[1;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(4,1);
dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
dC=dcdt;
end
C=Cv;
end
I estimate the initial conditions for the differential equations as the last ‘n’ values of the parameter vector, where ‘n’ are the number of differential equations in the systen, since it is easier to estimate them as parameters than to fix them as constants. It is also more robust.
I set ‘PopSz’ to 50 so it will run here in the time permitted. Increase it to at least 100 to get better initial results.
Use the ‘optshf’ options to use a hybrid function to fine-tune the estimated parameters using the fmincon function. (I use ‘optsAns’ here because some options I use offline are not permitted with the online Run feature.) Use ‘opts’ or ‘optshf’ instead offline.
.
6 个评论
Samhitha Tumkur
2022-10-22
Star Strider
2022-10-22
My pleasure!
‘for my "c" values, do i need a set of values for every equation of my dy function or are the last three good enough?’
You have to have one ‘dy’ equation for every column of ‘c’. that you want to fit. You can have more ‘dy’ equations of you want, providing they all share the same st of parameters, however the ones used to fit the data have to match the columns in ‘c’.
‘is there any way i can make this work so that i don't have to put a value for every equation in my 23 dy equations?’
I have no idea because I do not know what you are doing. As I mentioned, there must be one ‘dy’ equation for every column of ‘c’ that you want to fit. If you only want to fit three columns, you only need three differential equations. If all the equations in the system share the same parameters, and all those parameters can be estimated accurately by fitting three columns with three differential equations, then fitting three columns to estimate the parameters is good enough, providing they will accurately estimate all the necessary system parameters.
You will have to fit as many columns of ‘c’ as necessary to estimate all the necessary system parameters. However it would be best to fit all 23 of the colums with 23 separate ‘dy’ equations for the most optimal results.
.
Samhitha Tumkur
2022-10-31
Star Strider
2022-10-31
To return only the equations you want to fit, change this:
C=Cv;
so if you only wnated to return the last three (in this example), it might be:
C=Cv(:,[2 3 4]);
since the ordinary differential equation integration functions return column-wise solutions.
The error messages mean that the dimensions of the data and the objective function output do not match. If you only want to fit three columns of data, and your data have more columns than that, or you are fitting more columns than you have data for them, that error would be thrown.
Another possibility is that even if the column dimensions match, the row dimensions may not match if the differential equation integrating function encountered a singularity and stopped before it had finished the integration. In that event, just start the ga process over. If it keeps occurring, it may be necessary to look closely at the initial population to be certain that the parameter ranges are realistic.
I cannot tell from the posted error messages exactly what the problem is.
Samhitha Tumkur
2022-11-8
Star Strider
2022-11-8
As always, my pleasure!
更多回答（1 个）
Vinayak Choyyan
2022-10-14
0 个投票
Hi,
As per my understanding, you would like to implement a genetic algorithm to fit parameters.
Following are some resources that might help you with your issue.
- Genetic Algorithm - MATLAB & Simulink (mathworks.com)
- What Is a Genetic Algorithm? - Video - MATLAB (mathworks.com) - This contains how various approaches of genetic algorithm can be applied to given equations.
- Simple code for genetic algorithm - File Exchange - MATLAB Central (mathworks.com)- This contains a sample code which shows how genetic algorithm can be used to minimize or maximize equations.
类别
在 帮助中心 和 File Exchange 中查找有关 Genetic Algorithm 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
