Search code examples
matlabfminsearch

Error in fminsearch to solve this 4-variable problem?


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

Solution

  • 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