MATLAB: How to fit data on only one ordinary differential equations (out of multiple equations) with multiple unknown parameters

curve fitting

Hi,

I am trying to estimate 5 unknown parameters that exist in 10 ordinary differential equations (ODEs) by fitting it to experimental data of only one equation of the ODE.

I used the Monod Kinetics and curve fitting code to estimate the parameters but I got negative parameters instead of postive. This is the output.

Maximum number of function evaluations exceeded;

increase options.MaxFunEvals

Rate Constants:

beta(1) = 1.22325

beta(2) = -0.04436

beta(3) = 0.23296

beta(4) = -0.47970

beta(5) = -0.56506

This is my code

function S = Lassa(beta, t)

x0 = [100000;2000;400;400;5;50;200;500;50;20];

[T,Sv] = ode45(@DifEq, t, x0);

function dS = DifEq(t, x)

xdot = zeros(10,1);

piH = 0.000215; muH = 0.0000548;d1=0.01; piR=0.56;eR=0.85;rhoR=0.3;muR=0.6;d2=0.001;rho3=0.2;

rho6=0.0232;rho2=0.05;rho1=0.9;rho7=0.45;rho5=0.85;omega=0.45;rho4=0.06;gamma=0.03;

NR=x(8)+x(9)+x(10); NH=x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7);

xdot(1)=piH – beta(4)*x(9)*x(1)/NR -(beta(1)*x(5)+beta(2)*x(6)+beta(3)*x(7))*x(1)/NH -muH*x(1)+rho3*x(3)+rho6*x(7);

xdot(2)=beta(4)*x(9)*x(1)/NR+(beta(1)*x(5)+beta(2)*x(6)+beta(3)*x(7))*x(1)/NH -((1-omega)*rho2+ omega*rho1+muH)*x(2);

xdot(3)=(1-omega)*rho2*x(2)-(rho3*rho4+muH)*x(3);

xdot(4)=omega*rho1*x(2)-(rho5+muH)*x(4);

xdot(5)=rho5*x(4)-(rho7+d1+muH)*x(5);

xdot(6)=rho7*x(5)+rho4*x(3)-(rho6+d2+muH)*x(6);

xdot(7)=d1*x(5)+d2*x(6) – gamma*x(7);

xdot(8)=piR – beta(5)*x(9)*x(8)/NR – (muR+rhoR)*x(8);

xdot(9)=beta(5)*x(9)*x(8)/NR -(muR+rhoR+eR)*x(9);

xdot(10)=eR*x(9)-(muR+rhoR)*x(10);

dS = xdot;

end

S = Sv(:,5);

end

t=(1:37);

xdata = [1;5;2;5;7;4;1;0;2;9;3;1;0;1;0;0;0;0;1;0;0;0;0;0;0;0;0;0;3;0;0;0;0;0;0;0;1];

beta0=[0.5921;0.0162;0.01792;0.6506;0.7598];

[beta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@Lassa,beta0,t,xdata);

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

for k1 = 1:length(beta)

fprintf(1, '\t\tbeta(%d) = %8.5f\n', k1, beta(k1))

end

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

Cfit = Lassa(beta, tv);

Best Answer

The lsqcurvefit function allows you to set upper and lower bounds on the parameters it estimates with the ‘lb’ and ‘ub’ vectors. Set the lower bound to:

zeros(size(beta0))

and do not specify the upper bounds (set them equal to the empty array []) or set them equal to

Inf(size(beta0))

The original Monod Kinetics code did not need to constrain the parameters, so I did not constrain them.

Related Solutions

MATLAB: Errors when using the function ODE 45

tspan=[0 1];
% To run P-E system
%%Initialization
eta0=pi/4;
r0=sqrt(2);
beta0=0;
alpha0=pi/4;
xP0=0;
yP0=0;
xE0=1;
yE0=1;
ALPHA0=[eta0;r0;beta0;alpha0;xP0;yP0;xE0;yE0];
[tSol,ALPHASol]=ode45(@peFun,tspan,ALPHA0);      %function calling before function definition
plot(tSol,ALPHASol);
function fval=peFun(t,ALPHA)
% Get variables
eta=ALPHA(1);
r=ALPHA(2);
beta=ALPHA(3);
alpha=ALPHA(4);
xP=ALPHA(5);
yP=ALPHA(6);
xE=ALPHA(7);
yE=ALPHA(8);
% Constants
vER=3;
vPUR=5;
a=1;
% Describing dALPHA/dt
fval=zeros(8,1);
fval(1,1)=(-vPUR*sin(beta)+vER*sin(alpha))/r;
fval(2,1)=vER*cos(alpha)-vPUR*cos(beta);
fval(3,1)=(a/r)*(vER*sin(alpha)-vPUR*sin(beta));
fval(4,1)=(-vER*sin(alpha)+vER*sin(beta))/r;
fval(5,1)=vPUR*cos(eta+beta);
fval(6,1)=vPUR*sin(eta+beta);
fval(7,1)=vER*cos(eta+alpha);
fval(8,1)=vER*sin(eta+alpha);
end

MATLAB: How to estimate paramters.

I need help concerning parameter estimation and curve fitting using data. I am trying to estimate 22 paramters of the model using the the following matlab code. I keep geting this error in my work. Please I need help to correct my code. Thank you.

function lassafitw1(do_estimation)
warning off;
P_Data(:,1)=[1:80];
P_Data(:,2)=[ 0 16  117 252 305 386 454 525 574 646 717 783 848 910 1026 1261 1361 1424 1620 1673 1813 1940 2021 2304 2437  2950 3252 3410 3706 ...
    3896 5235 5338 4759 4862 5368 5586 6073 6190 6599 7109 7312 7897 8356 9004 9446 9780 10124 10934 11103 11466 11751 11841 11974 12022 12244 ...
    12313 12440  12593 127666 12735 12827 12884 12962 13012 13093 13150 13241 13402 13518 13586 13638 13701 13785 13894 13928 13945 13992 ...
    14001 14052 14122]./3800;
piH = 0.20;
muH = 0.02;
d1=0.2; 
piR=0.05;
eR=0.01;
rhoR=0.01;
muR=0.02;
d2=0.15;
rho3=0.05;
rho6=0.6;
rho2=0.6;
rho1=0.02;
rho7=0.2;
rho5=0.03;
omega=0.6;
rho4=0.5;
gamma=0.03;
beta1=0.2;
beta2=0.1;
beta3=0.08;
beta4=0.2;
beta5=0.02;
Sh0=100000; Eh0=2000; Eth0=0; Enh0=0; Ih0=0; Iqh0=0;Dh0=0; Sr0=500;Er0=0;Ir0=0;
INITIAL=[Sh0, Eh0, Eth0, Enh0, Ih0,Iqh0,Dh0, Sr0,Er0,Ir0];
Istart=1;    %month to start the model simulation
Iend=Istart + 81;
OPTIONS=odeset('AbsTol',0.001,'RelTol',0.001,'MaxStep',1/12);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Estimate parameters 
%by minimizing the sum of squares
%when fitting modeled to real prevalence data
    
do_estimation=1;
if(do_estimation)
xdata=P_Data(:,1)';
ydata=P_Data(:,2)';
x0(1,1)= 0.20;       %Lambda
x0(1,2)=0.02;      %mu
x0(1,3)= 0.2;  %beta
x0(1,4)= 0.05; %0.1;        %sigma
x0(1,5)=0.01; %0.005;        %gamma
x0(1,6)=0.01; %0.005;         %p
x0(1,7)= 0.02; %0.00001;        %theta
x0(1,8)= 0.15; %0.0008;        %delta1
x0(1,9)= 0.05; %0.004;         %delta2
x0(1,10)= 0.6;        %c
x0(1,11)= 0.6;     %omega
x0(1,12)= 0.02;     %alpha1
x0(1,13)=  0.2;     %alpha2
x0(1,14)=  0.03;     %alpha3
x0(1,15)=  0.6;     %alpha4
x0(1,16)= 0.5;   %alpha5
x0(1,17)= 0.03;
x0(1,18)= 0.2;
x0(1,19)= 0.1;
x0(1,20)= 0.08;
x0(1,21)= 0.2;
x0(1,22)= 0.02;
LB=[10 0 0.0001 5 0.0 0.000 0.0060 0.0060 0.0060 1 0.008  0.003 0 0 0 0 0 0 0 0 0 0 ];
UB=[50 1 1.5 20 0.9 0.9 0.9 0.9 0.9 10 1 1 1 1 1 1 1 1 1 1 1 5 5 ];
[x,Rsdnrm]=lsqcurvefit(@Model_Inc,x0,xdata,ydata,LB,UB,optimset);
'estimated parameters'
piH=x(1)
muH=x(2)
d1=x(3) 
piR=x(4)
eR=x(5)
rhoR=x(6)
muR=x(7)
d2=x(8)
rho3=x(9)
rho6=x(10)
rho2=x(11)
rho1=x(12)
rho7=x(13)
rho5=x(14)
omega=x(15)
rho4=x(16)
gamma=x(17)
beta1=x(18)
beta2=x(19)
beta3=x(20)
beta4=x(21)
beta5=x(22)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[t y] = ode45(@lassafitw1,[0:1/12:(Iend-Istart)], INITIAL);
Sh=y(:,1);
Eh=y(:,2);
Eth=y(:,3);
Enh=y(:,4);
Ih=y(:,5);
Iqh=y(:,6);
Dh=y(:,7);
Sr=y(:,8);
Er=y(:,9);
Ir=y(:,10);
incidence= beta1*y(:,5) + beta2*y(:,6);
%incidence= y(:,4)*gamma;

close all;
figure(1);
hold on
 h_l=plot(Istart+t,incidence,'r-');
set(h_l,'linewidth',2);
h_l=plot(P_Data(:,1),P_Data(:,2),'bo','Markersize',6);
set(h_l,'linewidth',2);
axis([Istart Iend 0 9]);
%ylim([0 45.00]) %                              The line added to change the y limit
%title('Ebola Prevalence(%)');
xlabel('Days','fontsize',15)
ylabel('Ebola Incidence','fontsize', 15)
hold off
function [ydot]=lassafitw1(t,y)
Sh=y(1);
Eh=y(2);
Eth=y(3);
Enh=y(4);
Ih=y(5);
Iqh=y(6);
Dh=y(7);
Sr=y(8);
Er=y(9);
Ir=y(10);
        ydot(1)=piH - beta4*y(9)*y(1)/(y(8)+y(9)+y(10)) -(beta1*y(5)+beta2*y(6)+beta3*y(7))*y(1)/(y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7))-muH*y(1)+rho3*y(3)+rho6*y(7);
        ydot(2)=beta4*y(9)*y(1)/(y(8)+y(9)+y(10))+(beta1*y(5)+beta2*y(6)+beta3*y(7))*y(1)/(y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7)) -((1-omega)*rho2+ omega*rho1+muH)*y(2);
        ydot(3)=(1-omega)*rho2*y(2)-(rho3+rho4+muH)*y(3);
        ydot(4)=omega*rho1*y(2)-(rho5+muH)*y(4);
        ydot(5)=rho5*y(4)-(rho7+d1+muH)*y(5);
        ydot(6)=rho7*y(5)+rho4*y(3)-(rho6+d2+muH)*y(6);
        ydot(7)=d1*y(5)+d2*y(6) - gamma*x(7);
        ydot(8)=piR - beta5*y(9)*y(8)/(y(8)+y(9)+y(10)) - (muR+rhoR)*y(8);
        ydot(9)=beta5*y(9)*y(8)/(y(8)+y(9)+y(10)) -(muR+rhoR+eR)*y(9);
        ydot(10)=eR*y(9)-(muR+rhoR)*y(10);
        ydot = ydot';
    
end % function 
function Inc=Model_Inc(x0,xdata)
% Inc=0;      %initialization of this not to have an empty array
piH=x0(1);
muH=x0(2);
d1=x0(3); 
piR=x0(4);
eR=x0(5);
rhoR=x0(6);
muR=x0(7);
d2=x0(8);
rho3=x0(9);
rho6=x0(10);
rho2=x0(11);
rho1=x0(12);
rho7=x0(13);
rho5=x0(14);
omega=x0(15);
rho4=x0(16);
gamma=x0(17);
beta1=x0(18);
beta2=x0(19);
beta3=x0(20);
beta4=x0(21);
beta5=x0(22);
%%% Initial values 
Sh0=100000; Eh0=2000; Eth0=0; Enh0=0; Ih0=0; Iqh0=0;Dh0=0; Sr0=500;Er0=0;Ir0=0;
INITIAL=[Sh0, Eh0, Eth0, Enh0, Ih0,Iqh0,Dh0,Sr0,Er0,Ir0];
% OPTIONS=odeset('AbsTol',0.0001,'RelTol',0.0001,'MaxStep',1/12);
% tdur=50;
[t y] = ode45(@lassafitw1, [0:1/12:Iend-Istart], INITIAL);
Sh=y(:,1);
Eh=y(:,2);
Eth=y(:,3);
Enh=y(:,4);
Ih=y(:,5);
Iqh=y(:,6);
Dh=y(:,7);
Sr=y(:,8);
Er=y(:,9);
Ir=y(:,10);   
incidence= beta1*y(:,5) + beta2*y(:,6)
%incidence= y(:,4)*gamma;
for i=1:length(xdata)
    
    ind=find(xdata(i)-Istart==t);
    Sh=y(ind,1);
    Eh=y(ind,2);
    Eth=y(ind,3);
    Enh=y(ind,4);
    Ih=y(ind,5);
    Iqh=y(ind,6);
    Dh=y(ind,7);
    Sr=y(ind,8);
    Er=y(ind,9);
    Ir=y(ind,10);
       
       Inc(i)= beta1*Ih + beta2*Iqh;  %incidence
      %Inc(i)= Is*gamma;  %incidence
 
end 
end     % end of Model_prev    
% Q1= gamma + mu; 
% Q2= mu + theta + delta1;
% Q3 = delta2 + sigma + mu;
% mstar= Lambda*gamma1 / omega* mu;
% R_M = beta*c*gamma*(rho + eta*sigma)*(1-mstar)*(p*theta + Q2*(1-p))/(rho*Q1*Q2*Q3)
end     %end of ebola_fit

Best Answer

Related Solutions

MATLAB: Errors when using the function ODE 45

MATLAB: How to estimate paramters.

Related Question