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

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

Answer 1

[Star Strider]

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.

[Chinwendu Madubueze]

Thank you Star Strider. I did that but i got the error message like this

??? Undefined function or variable 'beta0'.

Error in ==> Lassafit at 3

lb = zeros(size(beta0));

What should I do? How can I plot the data with the 5th equation of the ODE.

[Star Strider]

As always, my pleasure!

You posted this:

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

so I assume that you had already defined it.

Be sure to define ‘ydata’ as well as ‘xdata’, then put ‘lb’ here:

beta0=[0.5921;0.0162;0.01792;0.6506;0.7598];
lb = zeros(size(beta0));
[beta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@Lassa,beta0,t,xdata,ydata,lb,[]);

To return only the result of differential equation 5 as the output of the objective function, this:

S = Sv(:,5);

appears to be correct.

[Chinwendu Madubueze]

Thank you.

I think it is perfect except that some beta are zeros. This is the output

Optimization terminated: relative function value

changing by less than OPTIONS.TolFun.

Rate Constants:

beta(1) = 0.00000

beta(2) = 0.00000

beta(3) = 0.00000

beta(4) = 0.00000

beta(5) = 0.05030

Is there anything I can do so that the betas will not be zeros. Please I want to plot the data with the fifth equation in one plot.

Thank you for your help.

[Chinwendu Madubueze]

Please I want to plot the data with the fifth equation in one plot. How can I write the code for the plot.

Thank you for your help.

[Star Strider]

As always, my pleasure.

I have no idea what the probllem could be. Note that many of your ‘xdata’ values are zeros, (you did not post ‘ydata’) so it is certainly possible that the parameters are zeros. However, changing the format descriptor in the loop that prints those parameters:

for k1 = 1:length(beta)
    fprintf(1, '\t\tbeta(%d) = %23.15E\n', k1, beta(k1))
end

could give very small, non-zero results.

[Chinwendu Madubueze]

Thank you very much. I will try the code. You have save me from my code problem.

[Star Strider]

As always, my pleasure!