Genetic algorithm fot multi differential equations

good mornig i am working at a project in which i want to optimize the parameters of the 3 differential equations defined in the fun function, the parameters are 2, a = x(1) and b = x(2), i used the command ga to optimize them according to the available data that io already have, but i don't understand what i wrong because i get a lot of errors, also i don't know how to plot the final graph and compare it with the real one, here is the code:

clear all
N = 1e7;
dt = 1;
%% first  data curve
Tmax = 99;
t = 1:dt:Tmax;
data = xlsread('MyData.xlsx', 'Foglio1', 'B1:B100'); %infected people per day

%% Loading data into arrays
y = data(:,1); %copy data into array
I = y';  % vector of infected people per day
I0 = I(1);

%% setting conditions for the ga
A = []; 
b = [];
Aeq = [];
beq = [];
lb = [0 0];
ub = [5 5];
variables = 2;
%% defining the fitness function and the distance fitness function which will be inside the ga
fun = @(x) fitness_fun(x,Tmax,N,dt);
objFun = @(x) norm(fun(x)-I);

%% solution given by the ga 
coeff = ga(objFun,variables,A,b,Aeq,beq,lb,ub);

%% graph of the solution compared withj the real one
axes();
plot(t, I, 'b+');
hold on
plot(t, fun(coeff), 'r-');
legend({'Data points', 'Fitted Curve'})

the following is the fitness function which contains the ode45 command for the fun function with the differential equations:

function y = fitness_fun(x,Tmax,N,dt)

    tspan = 0:dt:Tmax;
    y0 = [N 1 0];
    [a,y] = ode45(@(x) fun, tspan, y0);

end

    %% series of function for the ODE solver
function dydt = fun(x,y)
    N = 1e4;
    a = x(1);
    b = x(2);
    dydt(1) = -a*y(1)*y(2)/N; %-a*S*I
    dydt(2) =  a*y(1)*y(2)/N - b*y(2); % a*S*I
    dydt(3) =  b*y(2);        % b*I
end

Answer 1

[Star Strider]

The ‘fitness_fun’ function should probably be something like this:

function y = fitness_fun(x,t,N)
    tspan = t;
    y0 = [N 1 0];
    [t,y] = ode45(@(x) fun, tspan, y0);

    %% series of function for the ODE solver
    function dydt = fun(t,y)
        dydt = zeros(3,1);
        N = 1e4;
        a = x(1);
        b = x(2);
        dydt(1) = -a*y(1)*y(2)/N; %-a*S*I
        dydt(2) =  a*y(1)*y(2)/N - b*y(2); % a*S*I
        dydt(3) =  b*y(2);        % b*I
    end
end

Also, the time (or independent variable) vector for the data must be the ‘tspan’ argument, or it will not be possible to compare the fitted result with the data. I added the zeros call so that ‘fun’ will return a column vector, and added ‘t’ as its first argument. It will inherit the parameter vector ‘x’ from the outer function workspace.

I cannot test this, however my changes should get it closer to working. (I generally pass the initial conditions vector as the last elements of the parameter vector so that the optimisation routine can fit them as well.)

Answer 2

[abdelghani msaad]

a the end your suggestion was usefull but i wasn doing the difference between the same dimension matrix

objFun = @(x) norm(I - fun(x));

I was a vector and fun(x) a matrix

[Star Strider]

They both have to have the same dimensions fo that to work.

[abdelghani takha]

hi sir

how to use Genetic algorithm for this differential equations

dx1dt = a0*x(1) - w0*x(2);

dx2dt = a0*x(2) + w0*x(1);

dx3dt = - sum(ai.*dti.*exp(-0.5*(dti./bi).^2))- 1.0*(x(3) - zbase);

[Star Strider]

Post this as a new Question.

Answer 3

[abdelghani msaad]

Here i prepare the array with the real data in might be something else in your case

    y = letturaExcel(t3(i),t3(i+1));
    %% defining the fitness function and the distance fitness function which will be inside the ga
    

Here i hold the fitness function in the variable fun, then i i hold in the variable objFun the real aim for my optimization which is the norm between the real data and the values calculated by the model in the fitness_fun

    fun = @(x) fitness_fun(x,t,N,y0);
    objFun = @(x) norm(fun(x) - y(:,2));
    

here you prepare some details to give to the ga so it will help it to work correclty form more info look at the ga documentation.

    %% solution given by the ga
    pop = 50;
    maxGen = 30;
    opts = optimoptions('ga', 'PopulationSize',pop, 'TolFun',1e-5,'MaxGenerations',maxGen, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
    

i save x wich is the array of coefficients optimized by ga and fval which is a value taht gives you an idea on how good is the optimization.

    [x,fval] = ga(objFun,variables,A,b,Aeq,beq,lb,ub,[],[],opts);
    [J,Jv] = fitness_fun(x,t,N,y0);
    % this last command you can decide...
    

here there is the fitness function

function [J,Jv] = fitness_fun(x,t,N,y0)

    [T,Jv] = ode45(@fun, t, y0);
    
    %differential equations to ode45
    function dY = fun(t,y)
        
        dydt = zeros(3,1);
        a = x(1);
        b = x(2);
        dydt(1) = -a*y(1)*y(2)/N; %-a*S*I/N
        dydt(2) =  a*y(1)*y(2)/N - b*y(2); % a*S*I/N - b*I
        dydt(3) =  b*y(2);        % b*I
        dY = dydt;
       
    end

J = Jv(:,2);

end

[abdelghani takha]

Thank you very much. I will give it a try and contact you if there are any problems.

[abdelghani takha]

@abdelghani msaad

hi sir

In the case of variables :

ai = [-60 -15 0 15 90];

bi = [1.2 -5 30 -7.5 0.75];

ci = [0.25 0.1 0.1 0.1 0.4];

How do I write Ib and Ub

[abdelghani msaad]

sorry for the late answer but could you post a new question where you write all the code that you have written so far? because it will be helpful for me also in order to help you in the right way

[abdelghani takha]

@abdelghani msaad

sorry for the late answer

%===My model===============================

xi = cos(ti);

yi = sin(ti);

ta = atan2(x(2),x(1));

dti = rem(ta - ti, 2*pi);

%%%My variables

ai = [1.2 -5 30 -7.5 0.75];

bi = [0.25 0.1 0.1 0.1 0.4];

ti = [-1.0472 -0.2618 0 0.2618 1.5708];

%%%My model

dx1dt = a0*x(1) - pi*x(2);

dx2dt = a0*x(2) + pi*x(1);

dx3dt = - sum(ai.*dti.*exp(-0.5*(dti./bi).^2))- 1.0*x(3) ;

Ib = [0 -0.25 0 0;-10 -0.1 -0.4 0;20 0 -1 0;-14 0 0 0;-1 0 0 0 ]; % lower limits for the variables

Ub = [2 0.4 0.5 1;10 0.2 0.4 1;40 0.2 0.1 1;0 0.20 0.4 1;1.5 0.8 3 1 ]; % upper limits for the variables

A = [];

b = [];

Aeq = [];

beq = [];

%==================================

nvar =3; % this is the number of variables in the objective functions

%==================================

% Now we call our objective function

[ECG,Y0] = fitness_fun(ti,ai,bi);

fun = @(x)fitness_fun(ti,ai,bi);

objFun = @(x) norm(fun(x) - I);

%fitnessF = @(ti,ai,bi)Mysim(ti,ai,bi);

%===================================

% Define your options here ,iterations,generations, etc

opt = optimoptions('ga', 'Display','iter',...% Display the iterations

'MaxGenerations', 1000*nvar, ...%Defins the number of generations% increasing improves accuracy

'PopulationSize', 50, ...%sets the populationsize

'FunctionTolerance', 1e-6, ...%this is the function tolerance % reducing improves accuracy

'PlotFcn', @gaplotbestindiv);% this command visualizes the results

%===========================================

%Now we optimise

% [x, fval] = ga(fitnessF,nvar,[],[],[],[],lb,Ub,[],opt);

% Note we have left some parts blank because we do not need them here

[x,fval] = ga(objFun,nvar,A,b,Aeq,beq,Ib,Ub,[],[],opt);

%% graph of the solution compared withj the real one

axes();

plot(t, I, 'b+');

hold on

plot(t, fun([x,fval]), 'r-');

legend({'Data points', 'Fitted Curve'})

[Omima Musa Mohmed Abusil]

I do have this code but the ga has failed at initialize the value and I do not know why? so will you help me to locate the issue and thanks

here is my code

function [J,Jv]=paramfun(theta,t,bt0)

% Monod Model for PPB growth

% dx/dt = Mumax*inhibition factor*x - kd*x

% ds/dt = -(1/y)*Mumax*inhibition factor*x

% with

% variable b(1) = x, b(2) = s

% parameter theta(1) = Mumax, theta(2) = inhibition factor, theta(3) = decay, theta(4)= yield

[T,Jv] = ode45(@fun,t,bt0);

function dC = fun(t,b);

dcdt = zeros(2,1);

dcdt(1)= theta(1)*theta(2)*b(1)-theta(3)*b(1);

dcdt(2) = -(1/theta(4))*theta(1)*theta(2)*b(1);

dC=dcdt;

end

J=Jv(:,2);

end

clear all

clc

t = [0 2.1 4 5 22.5 24.5 26.5 28.5 29.5 46.5 48.5 50.5 52.5 53.5 70.5 73.25 75.5]';

x = [19.5 23.57 24.33 24.33 80.6 100.76 142.2 174.14 188.21 321.3 331.18 324.33 322.432 322.432 327.755 332.052 319.77]';

s = [999.957 996.4 982.012 968.86 495.17 459.42 429.2 403.65 392.43 292.94 287.5 282.74 278.55 276.64 256.49 254.6 253.22]';

b = [x s];

% optimization

lb = [0,0,0,0];

ub = [0.04,0.5,0.03,1];

%%

A = [];

b = [];

Aeq = [];

beq = [];

%%

nvar = 4%;%number of variable in the objective function

%% to call the objective function

fun =@(b)paramfun(theta,t,bt0);

objFun = @(b)norm(fun(b)-b);

%fitnessF = @(theta,t,bt0)Mysim(theta,t,bt0);

%% defining the iteration option

options = optimoptions('ga','Display','iter','MaxGenerations',1000*nvar,'PopulationSize',55,'FunctionTolerance',1e-6,'PlotFcn',@gaplotbestindiv);

%% optimization

[b,fval]=ga(objFun,nvar,A,b,Aeq,beq,lb,ub,[],[],options);

%% drawing model data and real one

axes()

plot(t,b,'o'); hold on;

plot(t,fun([b,fval]),'-bo'); hold off;

waiting for your feedback urgently

[abdelghani msaad]

Hi man, im sorry but i am not an active user here, and it has been a while since i have not used more GAs, im not trained, what i can suggest you both is to post you questions as new Questions and not here in the comments so other people will be also able to see the and help you.

@Star Strider answers if the question is new to get prices in Mathworks so post a new question on the page.

Answer 4

[abdelghani takha]

@Star Strider

hi sir

In the case of variables :

ai = [-60 -15 0 15 90];

bi = [1.2 -5 30 -7.5 0.75];

ci = [0.25 0.1 0.1 0.1 0.4];

How do I write Ib and Ub

Answer 5

[abdelghani takha]

@Alan Weiss