MATLAB Answers

How to optimise this equation?

5 views (last 30 days)
Hey,
I have the following equation and I need to optimize this to find the best value for n:
Furthermore, I have the following data:
Tau = 5.2
E(t) = [0 0.0199867 0.0999334 0.1598934 0.1998668 0.1598934 0.1199201 0.0799467 0.05996 0.0439707 0.02998 0.011992 0]';
t = [0 1 2 3 4 5 6 7 8 9 10 12 14]';
Hopefully can someone help me to find the best way for optimizing this function, I have looked for some methods but didn't really knew how to apply them.
Thanks in advance,
Kind regards,
Danny
  3 Comments
Michal Kvasnicka
Michal Kvasnicka on 25 Sep 2019
So, do you mean by Optimization the best fit of the function E(t) on defined points, where n is optimized variable (1 <= n <= 20)? If yes, see answers. If no, please clarify.
As you can see, the optimized function in your case is not E(t), as you say, but
fun = @(n) norm(E - E_t(t,n,tau));

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 24 Sep 2019
Sinc you are using the factorial function, you are limited to integer values for ‘n’.
Try this:
objfcn = @(n,t,tau) (n.^n .* t.^(n-1)) ./ (factorial(n-1) .* tau.^n) .* exp(-n.*t/tau);
for k = 1:20
fcn(k) = norm(E - objfcn(k,t,Tau));
end
[Out,n] = min(fcn)
and the value of ‘n’ that produces the lowest residual norm ia 4.
If you want to use the gamma function instead:
objfcn = @(n,t) (n.^n .* t.^(n-1)) ./ (gamma(n) .* Tau.^n) .* exp(-n.*t/Tau);
[B,rsdnrm] = fminsearch(@(b) norm(E - objfcn(b,t)), 10)
evaluating to:
B =
4.2811
  2 Comments
Star Strider
Star Strider on 25 Sep 2019
My pleasure.
That is the initial parameter estimate for ‘n’. The fminsearch function (and all other nonlinear optimization functions that I am aware of) have to have a starting estimate for the parameters (here one parameter) the are required to estimate. I chose 10 here because it is in the centre of the range you specify.

Sign in to comment.

More Answers (1)

John D'Errico
John D'Errico on 24 Sep 2019
Edited: John D'Errico on 24 Sep 2019
Must n be an integer? After all, you have factorial(n-1) in there. If it must be an integer then just test a list of integer values for n. That will not be a long list, as you will have massive numerical problems if n is large. Pick the value of n that makes you happiest.
Hint: no, if n is an integer, you will never find an integer value for n that will come even remorely close to that data.
Instead, you very much need to use gamma(n) in there, instead of factorial(n-1). They are the same for integers, but the gamma function extends factorial nicely onto the real line.
Next, just to be safe, I used realpow in this. Although the .^ operator should be fine here.
Tau = 5.2
E = [0 0.0199867 0.0999334 0.1598934 0.1998668 0.1598934 0.1199201 0.0799467 0.05996 0.0439707 0.02998 0.011992 0]';
t = [0 1 2 3 4 5 6 7 8 9 10 12 14]';
E_t = @(t,n,tau) realpow(n,n).*realpow(t,n-1).*exp(-n*t/tau)./gamma(n)./realpow(tau,n);
Having done that, it is worth testing your fucntion to see if it has a chance to fit your data. So I played around with a few arbitrary values for n, and 4.25 seems to work resonably well, except that it does not reach the peak.
fplot(@(t) E_t(t,4.25,5.2),[0,14]),hold on,plot(t,E,'o')
As you see, it fits all ot the points pretty well, but it misses that peak. However, if you increase n to the point where is is large enough to hit the peak, then it looks like utter crap for the rest of the data.
fplot(@(t) E_t(t,6,5.2),[0,14]),hold on,plot(t,E,'o')
Could I have used an optimization to find the best value of n? Well, yes. That would be easy now. But the best value of n will be a compromise, likely close to 4.25, as I found initially.
fun = @(n) norm(E - E_t(t,n,tau));
[nopt,fval,exitflag] = fminsearch(fun,5)
nopt =
4.2811279296875
fval =
0.0301968802494492
exitflag =
1
fplot(@(t) E_t(t,nopt,5.2),[0,14]),hold on,plot(t,E,'o')
This points out several things. First, that I have a good eye for getting the best fit, even without a computer. But it also suggests that your data does not fit that model perfectly, missing that peak, even when the rest of the data fots nicely..

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by