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

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

[Star Strider]

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.

.

Answer 1

[Star Strider]

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

.

[Jaime Mata de la vega]

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

[Star Strider]

As always, my pleasure!