How to enter multiple parameter estimations in a single genetic algorithm (GA)

Question

0 个投票

sys_id2.m

Hello,

I have this code that estimates 2 variables (kd and Kla) from an ODE by using a genetic algorithm. The code predicts these two variables for a given data (time05 and C05). However, I have 3 data sets, 05, 20, and 35, and I want to find a kd and Kla that satisfy these 3 data sets.

There are 4 values in the ode. kd and KLa, which are the unknown values and that I want to find with the GA. O3_x which is a known value and is different for each data set as it is related to pressure. Finally, there is the initial condition, O3o which is always zero.

Does anyone know how I can integrate these 3 data sets into a single GA to get only one kd and one Kla that satisfy these 3 data sets?

% DATA

% set 1

time05 = [0; 50; 200; 300; 400; 600; 850; 1190; 1770; 2400];

C05 = [0;0.05e-04;0.25e-04;0.45e-04;0.55e-04;0.6e-04;0.75e-04;0.80e-04;0.89e-04;0.92e-04];

% set 2

time20 = [0;50;125;250;480;625;780;1000;1450;1700;2350];

C20 = [0e-03;0.0200e-03;0.0750e-03;0.1300e-03;0.1880e-03;0.2200e-03;0.2400e-03;0.2650e-03;0.3000e-03;0.3200e-03;0.3400e-03];

% set 3

time35 = [0;100;150;210;270;310;400;600;750;900;1220;1800;1910];

C35 = [0e-03;0.0600e-03;0.1300e-03;0.1600e-03;0.2200e-03;0.2500e-03;0.3000e-03;0.3400e-03;0.3750e-03;0.4250e-03;0.4600e-03;0.4850e-03;0.5050e-03];

% Optimoptions - Reproducibility

options = optimoptions('ga','ConstraintTolerance',1e-6,'PlotFcn', @gaplotbestf);

rng default % For reproducibility

% Genetic Algorithm

x = ga(@(par) sys_id2(par,time05,C05),3,[],[],[],[],[1e-4;0.0002;0],[0.0009;15;0],[],options);

Optimization terminated: average change in the fitness value less than options.FunctionTolerance.

% Parameters:

kla=x(1); kd=x(2);

% O3_x calculation by pressure

H = 5.22E+05; %

Pa05 = 500; % Pa

O3_x = Pa05/H;

Pa20 = 2000; % Pa

O3_x20 = Pa20/H;

Pa35 = 3500; % Pa

O3_x35 = Pa35/H;

% Initial condition

O3o=x(3);

% Definition of f(x,t)

fvdp = @(t,O3) -kd*O3^2+kla*(O3_x-O3);

% Solution

[t_05,y_05] = ode45(fvdp,time05, O3o);

% Plotting

plot(time05,C05,'-o','MarkerSize',6, ...

'MarkerEdgeColor','red', ...

'LineWidth',1, ...

'MarkerFaceColor','red', ...

'Color','red')

hold on

plot(t_05,y_05,'b--','LineWidth',2)

txt05 = {['KLa = ' num2str(x(1)) ' s-1'],['k = ' num2str(x(2)) ' M–1 s–1'],['O3* = ' num2str(O3_x) ' mol/L'],['P = ' num2str(Pa05) ' Pa']};

xt05 = [2000];

yt05 = [0.9e-4];

text(xt05,yt05,txt05)

set(gcf,'color','w');

title('Ozone absorption')

legend('Sotelo et al., 1989 0.5KPa','estimated 0.5 KPa')

xlabel('t (s)','Interpreter','Latex','FontSize', 12)

ylabel('$[O_3] (mol L-1)$','Interpreter','Latex','FontSize', 12)

hold off

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

Star Strider 2023-6-15

编辑：Star Strider 2023-6-16

在 MATLAB Online 中打开

I doubt that one set of kinetic parameters could estimate all three data sets, since they seem to be completely different —

% set 1

time05 = [0; 50; 200; 300; 400; 600; 850; 1190; 1770; 2400];

C05 = [0;0.05e-04;0.25e-04;0.45e-04;0.55e-04;0.6e-04;0.75e-04;0.80e-04;0.89e-04;0.92e-04];

% set 2

time20 = [0;50;125;250;480;625;780;1000;1450;1700;2350];

C20 = [0e-03;0.0200e-03;0.0750e-03;0.1300e-03;0.1880e-03;0.2200e-03;0.2400e-03;0.2650e-03;0.3000e-03;0.3200e-03;0.3400e-03];

% set 3

time35 = [0;100;150;210;270;310;400;600;750;900;1220;1800;1910];

C35 = [0e-03;0.0600e-03;0.1300e-03;0.1600e-03;0.2200e-03;0.2500e-03;0.3000e-03;0.3400e-03;0.3750e-03;0.4250e-03;0.4600e-03;0.4850e-03;0.5050e-03];

figure

plot(time05, C05, '.-', 'DisplayName','500 Pa')

hold on

plot(time20, C20, '.-', 'DisplayName','2000 Pa')

plot(time35, C35, '.-', 'DisplayName','3500 Pa')

hold off

grid

xlabel('Time')

ylabel('Concentration')

legend('Location','best')

title('Original Data')

C05i = interp1(time05, C05, time35, 'linear');

C20i = interp1(time20, C20, time35, 'linear');

t = time35;

Cmtx = [C05i C20i, C35];

figure

plot(t, Cmtx(:,1), '.-', 'DisplayName','500 Pa')

hold on

plot(t, Cmtx(:,2), '.-', 'DisplayName','2000 Pa')

plot(t, Cmtx(:,3), '.-', 'DisplayName','3500 Pa')

hold off

grid

xlabel('Time')

ylabel('Concentration')

legend('Location','best')

title('Interpolated Vectors')

The different pressures would need to be part of the kinetic equations (one for each pressure, I leave that to you), and they would all have to be interpolated to the same (shortest) time vector, with all original or interpolated concentrations at each of those times. Then, that would probably be possible. See: Estimate Parameters for System of ODEs with given Data using a Genetic Algorithm (COVID-19-Model) - MATLAB Answers - MATLAB Central for an example.

EDIT — (16 Jun 2023 at 17:42)

Added vector interpolations and plot.

.

请先登录，再进行评论。

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

Follow Question

Answer 1

Star Strider 2023-6-17

在 MATLAB Online 中打开

0 个投票

This works.

I will let you tweak it to give the appropriate results, since I have no idea what system you are modeling. I had to make several guesses as to what the code does and how it works (since it lacks significant documentation), so check to be certain I made the correct ones. I suggest increasing ‘PopSz’ to 250 or 500 to improve the possibility of getting the correct parametyer estimates. I kept it low here so the codee would run in tha allotted 55 seconds. It may require several different runs to get reasonable results.

% set 1

time05 = [0; 50; 200; 300; 400; 600; 850; 1190; 1770; 2400];

C05 = [0;0.05e-04;0.25e-04;0.45e-04;0.55e-04;0.6e-04;0.75e-04;0.80e-04;0.89e-04;0.92e-04];

% set 2

time20 = [0;50;125;250;480;625;780;1000;1450;1700;2350];

C20 = [0e-03;0.0200e-03;0.0750e-03;0.1300e-03;0.1880e-03;0.2200e-03;0.2400e-03;0.2650e-03;0.3000e-03;0.3200e-03;0.3400e-03];

% set 3

time35 = [0;100;150;210;270;310;400;600;750;900;1220;1800;1910];

C35 = [0e-03;0.0600e-03;0.1300e-03;0.1600e-03;0.2200e-03;0.2500e-03;0.3000e-03;0.3400e-03;0.3750e-03;0.4250e-03;0.4600e-03;0.4850e-03;0.5050e-03];

figure

plot(time05, C05, '.-', 'DisplayName','500 Pa')

hold on

plot(time20, C20, '.-', 'DisplayName','2000 Pa')

plot(time35, C35, '.-', 'DisplayName','3500 Pa')

hold off

grid

xlabel('Time')

ylabel('Concentration')

legend('Location','best')

title('Original Data')

C05i = interp1(time05, C05, time35, 'linear');

C20i = interp1(time20, C20, time35, 'linear');

t = time35;

Cmtx = [C05i C20i, C35];

figure

plot(t, Cmtx(:,1), '.-', 'DisplayName','500 Pa')

hold on

plot(t, Cmtx(:,2), '.-', 'DisplayName','2000 Pa')

plot(t, Cmtx(:,3), '.-', 'DisplayName','3500 Pa')

hold off

grid

xlabel('Time')

ylabel('Concentration')

legend('Location','best')

title('Interpolated Vectors')

ftns = @(parms) norm(Cmtx - kinetics(parms,t));

PopSz = 100;

Parms = 5;

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)

Start Time: 2023-06-17 20:03:01.4420

[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],optsAns);

Optimization terminated: average change in the fitness value less than options.FunctionTolerance.

t1 = clock;

fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)

Stop Time: 2023-06-17 20:03:54.6836

GA_Time = etime(t1,t0)

GA_Time = 53.2416

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)

Elapsed Time: 5.324160300000000E+01 00:00:53.241

fprintf('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)

Fitness value at convergence = 0.0006 Generations = 395

fprintf(1,'\tRate Constants:\n')

Rate Constants:

for k1 = 1:length(theta)

fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))

end

Theta( 1) = 0.00048 Theta( 2) = 35.88690 Theta( 3) = 0.00010 Theta( 4) = 0.00011 Theta( 5) = 0.00010

tv = linspace(min(t), max(t));

Cfit = kinetics(theta, tv);

c = Cmtx;

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(parms,t)

c0 = parms(3:5);

% Parameters:

kla = parms(1);

kd = parms(2);

% O3_x calculation by pressure

H = 5.22E+05; %

Pa05 = 500; % Pa

O3_x05 = Pa05/H;

Pa20 = 2000; % Pa

O3_x20 = Pa20/H;

Pa35 = 3500; % Pa

O3_x35 = Pa35/H;

[T,Cv] = ode15s(@(t,c)DifEq(t,c,kla,kd,O3_x05,O3_x20,O3_x35),t,c0);

% % Initial condition

% O3o=x(3);

% Definition of f(x,t)

function fvdp = DifEq(t,O3,kla,kd,O3_x05,O3_x20,O3_x35)

fvdp = zeros(3,1);

fvdp(1) = -kd*O3(1)^2+kla.*(O3_x05-O3(1));

fvdp(2) = -kd*O3(2)^2+kla.*(O3_x20-O3(2));

fvdp(3) = -kd*O3(3)^2+kla.*(O3_x35-O3(3));

end

C = Cv;

end

.

2 个评论
显示无隐藏无

Jaime Mata de la vega 2023-6-19

OMG @Star Strider Thank you very much. Yes, it is what I was looking for. I appreaciate it a lot

Star Strider 2023-6-19

As always, my pleasure!

请先登录，再进行评论。

How to enter multiple parameter estimations in a single genetic algorithm (GA)

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

采纳的回答

2 个评论
显示无隐藏无

更多回答（0 个）

类别

产品

版本

标签

Community Treasure Hunt

How to enter multiple parameter estimations in a single genetic algorithm (GA)

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

采纳的回答

2 个评论 显示 无 隐藏 无

更多回答（0 个）

类别

产品

版本

标签

另请参阅

Community Treasure Hunt

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

2 个评论
显示无隐藏无