How can I use levenberg-marquardt algorithm with two variables

Hello, everyone.

I started to use lsqcurvefit function in matlab for estimating the parameters in reaction kinetics.

My reaction kinetics are like this => r = k1[a][b] / (1+k2[a]+k3[b])^2

When I tried to use lsqcurvefit function, it demands only one variable x, however in my kinetics there are two variables a and b for output(reaction rate r).

Is there any way that can estimate reaction constant parameters k1,k2,k3 with two reaction variables a and b in kinetics?

I searched in matlab to solve this problem, but I couldn't find it.

I hope I could get good reply at here.

Thanks!

[J. Alex Lee]

what form is the data that you have (you have list of values for r with corresponding [a] and [b])? what does your model residual function look like, and how are you calling it?

based on what i see, you should be able to do it if you have more than 3 triplets of (r,a,b) data

[Taeksang Yoon]

Dear, J.Alex Lee

Thanks for answering my question.

1) In the reaction kinetics, a and b are the reactants of reaction, and [a] and [b] are the concentration of each reactants.

r is the reaction rates for the corresponding reactants concentration.

2) What does residual function mean? Is it necessary for estimation?

3)I have 12 data of (r,a,b) and I'm going to use 8 of them for parameter estimation and the rest will be used for validation. I think I need more data for accurate estimation.

Thanks!

[J. Alex Lee]

As @Star Strider's answer shows, I guess no, it is not important for you to know what the "residual function" is, because lsqcurvefit will automagically define it for you based on the model function that you give it. But it is definitely necessary for estimation since it is what defines "how far is your model from your data at this current set of estimates".

As for estimation and validation samples, I don't know if that's a typical thing for these types of goals, but if you're not going to use all available data I suppose I'd throw another caution out there to make sure you sample well...depending on how noisy your data is, 8 data points for 3 parameters doesn't seem like a lot, but i guess also depends on how "rigid" the model is...

[Taeksang Yoon]

Oh, I see.

I agree with your opinion. I think I need more experimental data, or modify my model with other assumptions.

Thanks for the reply!

Answer 1

[Star Strider]

It will be necessary to parameterize the ‘k’ constants (that you want to estimate) as:

k1 = p(1)
k2 = p(2)
k3 = p(3)

If the independent variable data are ‘a’ and ‘b’ and the dependent variable is ‘r’, then the data will be:

x = [a(:) b(:)];
y = r(:);

the objective function will be:

kinfcn = @(b,x) p(1).*x(:,1}.*x(:,2) ./ (1 + p(2).*x(:,1) + p(3).*x(:,2)).^2;

and the lsqcurvefit call would then be:

B = lsqcurvefit(kinfcn, rand(3,1), x, y, zeros(1,3))

The ‘zeros(1,3)’ prevents the kinetic constants (the ‘k’ parameters from being negative, since if I remember correctly, they should always be greater than 0.

[Taeksang Yoon]

Thanks for the detail instruction! It will help a lot!

I have one question in your answer.

In your code ' lsqcurvefit(kinfcn, rand(3,1),x,y,zeros(1,3))', does rand(3,1) indicates the random initial value of parameter k?

[Taeksang Yoon]

One more thing,

I tried to make a code as you instroduce, but when I run the code, the error comes out with the message 'The function value size and YDATA size are not the same.'

This is the code that I made with your guide and my matlab knowledge

clear all
clc

a=[0.68;0.51;0.34;1.02;0.77;0.51;1.36;1.02];
b=[0.68;0.68;0.68;1.02;1.02;1.02;1.36;1.36];
r=[307.2;283.5;239.0;435.6;407.5;341.2;542.8;499.5];

x=[a(:) b(:)];
y=r(:);



B=lsqcurvefit(@kinfcn, rand(3,1),x,y,zeros(1,3))

k1=p
k2=q
k3=r

and the code of function kinfcn is like this

function kinfcn=f(p,x)
kinfcn= @(b,x) p*x(:,1).*x(:,2)./(1+q*x(:,1)+r*x(:,2)).^2;
end

Do you know why the matlab returns the error?

[Star Strider]

There was an error in ‘kinfcn’ that I did not catch, however it was not the one you used. The one you used does not work.

With those errors corrected:

a=[0.68;0.51;0.34;1.02;0.77;0.51;1.36;1.02];
b=[0.68;0.68;0.68;1.02;1.02;1.02;1.36;1.36];
r=[307.2;283.5;239.0;435.6;407.5;341.2;542.8;499.5];

x=[a(:) b(:)];
y=r(:);

kinfcn = @(p,x) p(1).*x(:,1).*x(:,2) ./ (1 + p(2).*x(:,1) + p(3).*x(:,2)).^2;

B=lsqcurvefit(kinfcn, rand(3,1),x,y,zeros(1,3))

figure
plot3(a, b, r, 'p')
hold on
plot3(a, b, kinfcn(B,x), '-r')
hold off
grid on
xlabel('a')
ylabel('b')
zlabel('r')
title(sprintf('$r = \\frac{%.3f[a][b]}{(1+%.3f[a]+%.3f[b])^2}$', B), 'Interpreter','latex', 'FontSize',18)

producing:

The fit of the funciton to the data is quite good.

[Taeksang Yoon]

Oh, I got it.

Thanks for the corrected code and the result!

I also get the estimation result by using curve fitting toolbox, and it also shows the almost same result.

[Star Strider]

As always, my pleasure!

I do not have the Curve Fitting Toolbox, since I can do everything I need to do with the Statistics and Machine Learning Toolbox and the Optimization Toolbox and Global Optimization Toolbox.