fitting "K" paremeter in a adsorption breakthrough curve data

I haves this code

% This code is an example of the use of ODE in matlab. In this example they use a simple formula and a FDA method

% and the frenluich isoterm.

clc

clear

close all

%input data

L=0.2;

eps=0.8;

u=0.016;

kf=6e-06;

c0=0.5;

rho=980;

De=1.2e-09;

ap=2e-3;

n=1.31;

K=0.7; %adjust parameter

%length

Nz=51;

z=linspace(0,L,Nz);

dz=z(2)-z(1);

%time

t=[0:125:5000];

%initial condition

ICA=zeros(1,Nz); %PDE

ICB=zeros(1,Nz); %ODE

IC = [ICA,ICB];

%ODE solver

[t y]=ode15s(@f,t,IC,[],eps,u,kf,c0,rho,Nz,dz,De,ap,K,n);

%Define value

c=y(:,1:Nz);

q=y(:,Nz+1:2*Nz);

%Recalculation

c(:,1)=c0;

c(:,end)=(4*c(:,end-1)-c(:,end-2))./3;

%visual

cp=c(:,end)./c0;

plot(t,cp,'or')

grid on

title('Breakthrough curve')

xlabel('time (min)')

ylabel('c/c0')

axis([0 5000 0 1.2])

hold on

%function

function dydt=f(t,y,eps,u,kf,c0,rho,Nz,dz,De,ap,K,n)

dydt=zeros(length(y),1);

dcdt=zeros(Nz,1);

dqdt=zeros(Nz,1);

%Define value

c=y(1:Nz);

q=y(Nz+1:2*Nz);

%Boundary conditions

c(1)=c0;

c(end)=(4*c(end-1)-c(end-2))./3;

%interior

for i=2: Nz-1

qstar(i) =K.*c(i).^n;

dqdt(i)=(3.*kf./ap).*(qstar(i)-q(i));

dcdz(i)=(c(i+1)-c(i-1))./2./dz;

dc2dz2(i)=(c(i+1)-2*c(i)+c(i-1))./dz.^2;

dcdt(i)= De.*dc2dz2(i)-u*dcdz(i)-rho*((1-eps)./eps).*dqdt(i);

end

dydt=[dcdt;dqdt];

end

What I want to do is to be able to adjust the value of K to the experimental data I add. My intention is to add experimental time and concentration data. With this data adjust the K value to match the curves.

But I don't know how to do it, I don't know how to integrate the fminsearch in this code.

Can you help me please??

Answer 1

[Star Strider]

Estimating ‘K’ is likely straightforward, and requires creating an objective function similar to the ‘kinetics’ function here, defining ‘theta’ as ‘K’ followed by a vector of the initial conditions. The objective function integrates the differential equations, and returns an array the same size as the data that lsqcurvefit (in my code here) uses the estimate ‘K’.

In any event, having the data is necessary to provide a specific solution.

I would do something like this, although creating the independent (‘t’) and dependent (‘c’) variable arrays from the data read from the file —

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]; 

theta0 = rand(10,1);

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)));

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

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

	Rate Constants:

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

		Theta( 1) =  0.76483
		Theta( 2) =  0.23492
		Theta( 3) =  0.20878
		Theta( 4) =  0.49180
		Theta( 5) =  0.62211
		Theta( 6) =  0.00000
		Theta( 7) =  0.90287
		Theta( 8) =  0.07146
		Theta( 9) =  0.02840
		Theta(10) =  0.00000

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','best')

function C=kinetics(theta,t)
    c0=theta(end-3:end);
    [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

.