How do I solve two ode (dct/dt and dR/dt) simultaneously with ODE45 and optimize two parameters using global optimization ?

2 views (last 30 days)
Gokul Rao
Gokul Rao on 7 Jun 2021
function dYdt = concentration(t,Y,p)
%Known Parameters
r= 250; %0.684 um radius of adsorbent
rf= 0.684; %Radius of concentration front of solute penetrating adsorbent
P=2.439; %2.439 density of adsorbent
m=1; %1 g mass of adsorbent
g=50; %50 mL/min gas flowrate
CO=0.1; %0.1 initial CO2 concentration
qet=138.72; %138.72 CO2 concentration at equilibrium
kF=43.764; %43.764 Freundlich constant
%Unknown Parameters
kg = p(1);
De = p(2);
%nonlinear differential equation to show change in concentration vs time
Ch= (qet*m)/(CO*g);
Bi= (kg/De)*r;
CET= cet/CO;
R= rf/r;
N= (3*Ch*qet*r^2)+((kF/nF)*(CET^(1/nF-1))*(Ch*Bi*(1-r^3)*CT)*(1/(r+Bi*(1-r))^2));
M= 1+Ch*(1-r^3)*(kF/nF)*(CET^(1/nF-1))*(Bi*(1-r))*(1/(r+Bi*(1-r)));
dRdt= -(Bi*(CO/P*qet)*(CT-CET))*(1/R^2);
dCTdt= N/M*(dRdt);
dYdt= [dCTdt;dRdt];
% Define the objective
function obj = objective(p)
%Initial condition at t=0
Y0= [R;CT];
%Specific time points corresponding to experimental time points (minutes)
tspan= [0 5 10 15 20 25 30 35 40 45 50 55 60];
%function handle to pass p through substrate function and obtain numerical
%model by solving nonlinear differential equation
[t,Y]= ode45(@(t,Y)concentration(t,Y,p),tspan,Y0);
%Experimental data at each time points ts
CT1measured= [1 2 3 4 5 6 7 8 9 10 11 12 13];
CT1measured= CT1measured';
%Objective function to minimise the square of the difference between the
%numerical model and experimental data
A = (CT-CT1measured).^2;
obj = sum(A);
close all;
%Parameter Initial Guess
p0= [kg;De];
%Objective function is the square of the difference between numerical and
fun = @objective; % function handle created
%identify unknown parameters
%[p,fval], find both the location of the minimum p of the nonlinear
%function and the value of the function at that minimum.
[p,fval] = fminunc(fun,p0);
%Optimised parameter values
kg = p(1);
De = p(2);
disp(['kg: ' num2str(kg)])
disp(['De: ' num2str(De)])
%Calculate model with updated parameters
p = [kg;De];
tspan= [0 5 10 15 20 25 30 35 40 45 50 55 60];
Y0= [R;CT];
[t,Y]= ode45(@(t,Y)concentration(t,Y,p),tspan,Y0);
CT1measured= [1 2 3 4 5 6 7 8 9 10 11 12 13];
%Plot of updated parameter values
plot(tspan,CT1measured,'o',tspan,Y,'r');xlabel('time (min)');ylabel('Ct/Co');
title('Ct/Co vs time (min)');grid on

Sign in to comment.

Answers (3)

Sulaymon Eshkabilov
Sulaymon Eshkabilov on 5 Jun 2021
These are coupled ODEs. To solve these ODEs numerically:
(0) Specify all constant parameters
(1) Build fcn file or anonymous fcn handle
(2) Call and solve the fcn using one of these solvers: ode23, ode45, ode113

Alan Weiss
Alan Weiss on 6 Jun 2021
For a relevant example, see Optimize an ODE in Parallel.
Alan Weiss
MATLAB mathematical toolbox documentation

Star Strider
Star Strider on 7 Jun 2021

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by