0

I am trying to solve the equation below using fminsearch, but I think the objective function is wrong.

How should I write the objective function or modify any other part of the code? This is basically a fitting problem, where the optimization procedure should fit the equation to the given data.

% Consider the following data:

Data = ...
  [0.0000    5.8955
   0.1000    3.5639
   0.2000    2.5173
   0.3000    1.9790
   0.4000    1.8990
   0.5000    1.3938
   0.6000    1.1359
   0.7000    1.0096
   0.8000    1.0343
   0.9000    0.8435
   1.0000    0.6856
   1.1000    0.6100
   1.2000    0.5392
   1.3000    0.3946
   1.4000    0.3903
   1.5000    0.5474
   1.6000    0.3459
   1.7000    0.1370
   1.8000    0.2211
   1.9000    0.1704
   2.0000    0.2636];

% Let's plot these data points.
t = Data(:,1);
y = Data(:,2);

plot(t,y,'ro')
title('Data points')
hold on

% fit the function: y =  c(1)*exp(-lam(1)*t) + c(2)*exp(-lam(2)*t)
%
% define the parameters in terms of one variable x:
%  x(1) = c(1)
%  x(2) = lam(1)
%  x(3) = c(2)
%  x(4) = lam(2)
%
% Then define the curve as a function of the parameters x and the data t:

F = @(x,t)(x(1)*exp(-x(2)*t) + x(3)*exp(-x(4)*t));

% We arbitrarily set our initial point x0 as follows: c(1) = 1,
% lam(1) = 1, c(2) = 1, lam(2) = 0:

x0 = [1 1 1 0];

% We run the solver and plot the resulting fit
options = optimset('TolFun',1e-5,'TolX',1e-5,'MaxFunEvals',10,'MaxIter',4000,'Display','iter');
[x,fval,exitflag,output] = fminsearch(F,x0,options)

plot(t,F(x,t))
hold off
ShadowWarrior
  • 180
  • 1
  • 12

1 Answers1

2

You're right, your objective function doesn't make sense. You could do a least squares fitting. Then the objective function should be defined as:

F = @(x,t) (x(1)*exp(-x(2)*t) + x(3)*exp(-x(4)*t));
Obj = @(x) (norm(F(x, Data(:,1))-Data(:,2)));

Then

x0 = [1 1 1 0];
options = optimset('TolFun',1e-5, 'TolX', 1e-5, 'MaxFunEvals',1000, 'MaxIter', 4000,'Display','iter');
[x,fval,exitflag,output] = fminsearch(Obj,x0,options);

tp = 0:0.01:2;
plot(Data(:,1), Data(:,2), 'ro');
title('Data points')
hold on
plot(tp,F(x,tp));
hold off

gives me:

enter image description here

EDIT:

Assuming you already know a parameter and you want to use function handles, you could do it like this

p = ... % Your calculation to get the parameter. In your case x(3) from the previous F
F = @(x, p, t) (x(1)*exp(-x(2)*t) + p*exp(-x(3)*t));
helper = @(x, p) (norm(F(x, p, Data(:,1))-Data(:,2)));
Obj = @(x) (helper(x, p));

x0 = [1 1 0]; % Note that there's now one variable/parameter less
options = optimset('TolFun',1e-5, 'TolX', 1e-5, 'MaxFunEvals',1000, 'MaxIter', 4000,'Display','iter');
[x,fval,exitflag,output] = fminsearch(Obj,x0,options);

tp = 0:0.01:2;
plot(Data(:,1), Data(:,2), 'ro');
title('Data points')
hold on
plot(tp,F(x,p,tp)); % Note that you need to pass p to F
hold off
joni
  • 6,840
  • 2
  • 13
  • 20
  • Thank you, is there some more accurate/robust fitting possible beside least squares? – ShadowWarrior Sep 30 '18 at 03:04
  • If you have the curve fitting toolbox installed, you could take a look at [fit](https://www.mathworks.com/examples/curvefitting/mw/curvefit-ex56454479-robust-fitting) which has a option for robust fitting in the sense that outliers will be lower weighted in the fit. But for your example data there aren't any outliers, so the least squares fitting should be fine. On the other side you could do a [spline interpolation](https://de.mathworks.com/help/matlab/ref/spline.html). – joni Sep 30 '18 at 06:25
  • there is also a build in particle swarm optimization for matlab – albusSimba Sep 30 '18 at 06:36
  • @JohnnyDrama, Thank you. Can you please suggest how F and obj would change when x(3) will come from ode45? – ShadowWarrior Sep 30 '18 at 17:46
  • @JohnnyDrama, this has been immensely helpful, thank you! I tried to adapt the above code to a new problem where parameters inside ode45 need to be fit and I could not make it work. Would it be possible for you to kindly have a look at the code? I have pasted it at codeshare and it can be edited there - https://codeshare.io/29Z4Y4 – ShadowWarrior Oct 01 '18 at 19:45
  • @ShadowWarrior Took a look at it and there are a few things that are unclear to me. I would recommend to post a new question here and provide a minimal example that exactly describes what you're trying to do. – joni Oct 01 '18 at 20:23
  • @JohnnyDrama, I have solved the problem. Just one question, In `(norm(F(x, p, Data(:,1))-Data(:,2)).*2)` why we are multiplying it by 2? Should not `(norm(F(x, p, Data(:,1))-Data(:,2))) = sqrt(sum(F(x, p, Data(:,1))-Data(:,2)).^2) = least squares method`?? – ShadowWarrior Oct 14 '18 at 06:38
  • You're right, it's not necessary to multiply it by factor 2, it's just a typo. Thanks for pointing out. :) But it won't change the result since it's just a constant factor. – joni Oct 14 '18 at 06:50
  • @JohnnyDrama, thank you for helping me out, highly appreciate it. – ShadowWarrior Oct 14 '18 at 09:15