Parameter estimation in ODE

I tried to use Star Strider ‘Igor_Moura’ function to solve same problem adding another one ode equation (without experimental data, just need to solve), but i have problem:

function C=kinetics(theta,t)
c0=[1;0;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(5,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);
    dcdt(5)=theta(1).*c(2)-theta(5);
    dC=dcdt;
    end
C=Cv;
end

So, I just added 5th equation, for example:

dcdt(5)=theta(1).*c(2)-theta(5)

I assuemd that there is no experimental data for 5th equation, and i just have to solve this equation.

I change nothing in script.

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=[1;1;1;1;1;1];

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

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

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'Location','N')

And i get this:

Error using lsqcurvefit (line 251)

Function value and YDATA sizes are not equal.

Error in script (line 28)

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

I would really appreciate the help.

Thanks in advance,

Answer 1

[Star Strider]

You are trying to fit 5 columns of your differential equation to 4 columns of data. That will throw the error you got.

Change the ‘C’ assignment in the ‘kinetics’ objective function to:

C=Cv(:,1:4);

and it works.

So ‘klinetics’ is now:

function C=kinetics(theta,t)
c0=[1;0;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(5,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);
    dcdt(5)=theta(1).*c(2)-theta(5);
    dC=dcdt;
    end
C=Cv(:,1:4);
end

If you then want to estimate Compartment 5, run your differential equation function again with all the estimated parameters to include it and plot it.

I created a version of ‘kinetics’ called ‘kinetics5’ to calculate and plot Compartment 5 as well. That entire additional code (with all the previously estimated parameters, so all this goes after the lsqcurvefit call) is:

function C=kinetics5(theta,t)
c0=[1;0;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(5,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);
    dcdt(5)=theta(1).*c(2)-theta(5);
    dC=dcdt;
    end
C=Cv;
end

for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics5(theta, tv);
figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'C_5(t)', 'Location','SW')

That should do what you want.

[Bosnian Kingdom]

Thank you for the answer:

I do this:

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=[1;1;1;1;1;1];

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

function C=kinetics5(theta,t)
c0=[1;0;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(5,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);
    dcdt(5)=theta(1).*c(2)-theta(5);
    dC=dcdt;
    end
C=Cv;
end
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics5(theta, tv);
figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'C_5(t)', 'Location','SW')

and gettng this error:

Error: File: script.m Line: 45 Column: 1

Function definitions in a script must appear at the end of the file.

Move all statements after the "kinetics5" function definition to before the function definition.

What is the problem?

[Bosnian Kingdom]

I do this, I made two separate function ypu write to me:

function C=kinetics(theta,t)
c0=[1;0;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(5,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);
    dcdt(5)=theta(1).*c(2)-theta(5);
    dC=dcdt;
    end
C=Cv(:,1:4);
end

and

function C=kinetics5(theta,t)
c0=[1;0;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(5,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);
    dcdt(5)=theta(1).*c(2)-theta(5);
    dC=dcdt;
    end
C=Cv;
end

and than i made script

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=[1;1;1;1;1;1];

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

tv = linspace(min(t), max(t));
Cfit = kinetics5(theta, tv);
figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'C_5(t)', 'Location','SW')

and got the results, but there is no writing thetas in command windows and I just get plot. Walues of thetas must found in Workspace, how to write me results for thetas in Command Windows however, i also got this:

Local minimum found.

Optimization completed because the size of the gradient is less than

the default value of the optimality tolerance.

<stopping criteria details>

[Bosnian Kingdom]

And also i cant' se results for the dcdt(5).

[Star Strider]

Save ‘kinetics’ as ‘kinetics.m’ and ‘kinetics5’ as ‘kinetics5.m’ on your MATLAB user path.

See the documentation on Function Basics (link).

[Bosnian Kingdom]

Thank you i solved th problem, however i still getting:

Local minimum found.

Optimization completed because the size of the gradient is less than

the default value of the optimality tolerance.

<stopping criteria details>

And i got all the results and plot.

[Star Strider]

As always, my pleasure.

The notice simply means that the solver found a minimum, although based on the value of the function at that point, it may not be the global minimum. The parameter set I got was:

		Theta(1) =  0.89197
		Theta(2) =  0.26069
		Theta(3) =  0.24870
		Theta(4) =  0.48183
		Theta(5) =  0.60977
		Theta(6) = -0.01301

You can also constrain the parameter estimates using lsqcurvefit so that none of them are negative. That is likely appropriate for a kinetic model.

[Bosnian Kingdom]

Thank you for your answer, can you explain me how to constrain the parameter estimates using lsqcurvefit so that none of them are negative. I solving kinetic model similar to this.

Is it posible in this problem you help me to solve (optimize kinetic parameter, solve 5th rquation set objective (for example that dcdt(1)/dcdt(2)=2) and then to solve problem?

function dC=DifEq(t,c)
    dcdt=zeros(5,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);
    dcdt(5)=theta(1).*c(2)-theta(5);
    dC=dcdt;

Thank in advance?

[Star Strider]

Add the ‘lb’ (lower bound) and ‘ub’ (upper bound) parameter constraints:

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

‘Is it posible in this problem you help me to solve (optimize kinetic parameter, solve 5th rquation set objective (for example that dcdt(1)/dcdt(2)=2) and then to solve problem?’

I cannot imagine how to do that, other than replacing ‘dcdt(2)’ with:

dcdt(2) = dcdt(1)*2;

That would likely work (I did not test it), although it may cause serious problems with the parameter estimation, unless that is also true for the data you are fitting. You would most likely need to change your model to reflect that.

[Bosnian Kingdom]

Thanks, it is working when i write this:

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

is it [0.1] lower band?

[Star Strider]

As always, my pleasure.

‘is it [0.1] lower band?’

No. That is the uppper bound, likely of only the first parameter. You do not need upper bounds for your parameters. (They will likely all be less than 1, however I do not bound them so that I can be certain my code is working correctly. If any parameters are greater than 1, I search for the reason, since it likely indicates an error in my code.) I included an empty matrix [] in the argument list for the upper bound simply for illustration (and since I always supply arguments for both bounds in my code). Either use the empty matrix or do not supply the upper bound argument at all, if you do not want to set an upper bound.

[Bosnian Kingdom]

Thank you, you are very helpful!

You've solved my dilemmas about this example

[Star Strider]

As always, my pleasure.

Thank you.