Search code examples
matlabconstraintscurve-fittingfminsearch

matlab: constrained fit with lsqcurvefit and fmincon


I have a set of data points (data_x, data_y). I need to fit a model function into this data. Model is a function of 5 parameters, and I have defined it like that:

function F = model(x,xdata)

fraction1 = x(4);
fraction2 = x(5);
fraction3 = 1-x(4)-x(5);

F=1-(fraction1.*(exp(-(xdata)./x(1)))+(fraction2.*(exp(-(xdata)./x(2))))+(fraction3.*(exp(-(xdata)./x(3)))));

parameters x(4) and x(5) are used to define three fractions, so their sum MUST be 1. To fit this function I was using lsqcurvefit, like that:

%% initial conditions
a0 = [guess1 guess2 guess3 0.3 0.3];

%% bounds
lb = [0 0 0 0 0 ];
ub = [inf inf inf 1 1];

%% Fitting options
curvefitoptions = optimset( 'Display', 'iter' );

%% Fit
a = lsqcurvefit(@model,a0,x,y,lb,ub,curvefitoptions);

The thing is that don't know how to add constraints, to keep the sum of fractions = 1. I know that lsqcurvefit is not the best solution for this problem, but I have no idea how to feed fmincon with these data to find my parameters. Many thanks for help!

EDIT: just a note, the max value of F could be 1... I was trying to cheat i.e. by adding or even multiplying everything by something like 10^(1-fraction1-fraction2-fraction3), but then I was ending up with almost equal fractions (0.33), what has no sense at all, cause other parameters are screwed... When I was fitting the same data using Origin (with same model + constrains) it works perfectly... When I used fixed Origin output fraction parameters fits were also great, but... its not the way to do that having dozen of fits to do :(


Solution

  • If you use fmincon for this (and use another parameter as the third fraction) the constraints are quite simple. You may need to play around with the fmincon options to get good convergence.

    function solution = my_fit_fun(xdata, ydata, a0)
    
    lb = [0 0 0 0 0 0];
    ub = [inf inf inf 1 1 1];
    
    %Aeq and beq specify that the last three parameters add to 1
    Aeq = [0 0 0 1 1 1];
    beq = 1;
    
    solution = fmincon(@objective,a0,[],[],Aeq,beq,lb,ub);
    
        function F = model(x)
    
            fraction1 = x(4);
            fraction2 = x(5);
            fraction3 = x(6);
    
            F=1-(fraction1.*(exp(-(xdata)./x(1)))+(fraction2.*(exp(-(xdata)./x(2))))+(fraction3.*(exp(-(xdata)./x(3)))));
    
        end
    
        function f = objective(x)
    
            yfit = model(x);
            f = sum((yfit - ydata) .^2);
        end
    
    end