How to constraint the values of fitted parameters with lsqcurvefit?

I am fitting a multiexponential function to a data set, but the solutions that lsqcurvefit is finding are away from realistic values. If I wanted to constraint the value of the parameters, to be greater than 0 but to not be grater that other value, lets say 0.001, how would I do it?

I will include my code below:

clear; clc; clf; close all;

D_0 = [0.5 0.00001 0.5 0.00001]; %initial guess
xdata = [9.85	32.05	66.83	114.44	174.55	247.57	333.00	431.44	542.19	666.04	802.12	951.39	1113.38	1287.47	1474.87	1674.29	1887.11	2600.14	2863.78	3139.17	3428.23	4043.41	4369.45	4709.33	5425.99	5802.67	6193.38	6595.39];
ydata = [1	0.9569	0.9528	0.8894	0.8387	0.8995	0.7911	0.773	0.7523	0.7155	0.7086	0.6478	0.6269	0.6175	0.574	0.551	0.4991	0.4559	0.4449	0.4314	0.4212	0.407	0.3856	0.3511	0.3526	0.303	0.3148	0.2912];

fun = @(D,xdata) D(1)*exp(-xdata.*D(2))+ D(3)*(exp(-xdata.*D(4)));
D = lsqcurvefit(fun, D_0, xdata, ydata);

semilogy(xdata, ydata, 'ko', xdata, fun(D,xdata), 'b-')
legend('Data', 'Fit')

format short
D_1 = D(1)
D_2 = D(2)
D_3 = D(3)
D_4 = D(4)

% I need to constraint the 2nd and 4th parameter to be below 0.001, even if
% it is not the most perfect fit.

% I this case D_2 is greater than 0.001

Answer 1

[Alex Sha]

hi, the result is good enough

Sum Squared Error (SSE): 0.0105245967805521
Root of Mean Square Error (RMSE): 0.01938758511131
Correlation Coef. (R): 0.99609231161877
R-Square: 0.992199893266025

Parameter	Best Estimate
----------	-------------
d1	0.364389939850661
d2	0.00133809356253748
d3	0.610149012103863
d4	0.000110756764815313

[Alfredo Scigliani]

I know that the fitting is good, but it is a false minimum. This fitting tells me about properties of a material in which something above 0.001 is just non sensical. There shold be another solution to the data fitting that the function is not getting to because it is stopping the iterations at this one.

Answer 2

[Mathieu NOE]

hello

I admit, my solution is the "poor man" solution as I don't have the optimization toolbox. But even with fminsearch I could do a 3 parameters fit that looks ok (fig 2, right ) by forcing the last parameter (d) to be = 0.001. If you let the 4 parameters free, d would slightly overshoot the limit (0.0013) , as displayed in fig 1 .

visualy , both solutions seems ok

clear; clc; clf; close all;
D_0 = [0.5 0.00001 0.5 0.00001]; %initial guess
xdata = [9.85	32.05	66.83	114.44	174.55	247.57	333.00	431.44	542.19	666.04	802.12	951.39	1113.38	1287.47	1474.87	1674.29	1887.11	2600.14	2863.78	3139.17	3428.23	4043.41	4369.45	4709.33	5425.99	5802.67	6193.38	6595.39];
ydata = [1	0.9569	0.9528	0.8894	0.8387	0.8995	0.7911	0.773	0.7523	0.7155	0.7086	0.6478	0.6269	0.6175	0.574	0.551	0.4991	0.4559	0.4449	0.4314	0.4212	0.407	0.3856	0.3511	0.3526	0.303	0.3148	0.2912];

%% 4 parameters fminsearch optimization
f = @(a,b,c,d,x) a.*exp(-b.*x) + c.*exp(-d.*x);
obj_fun = @(params) norm(f(params(1), params(2), params(3), params(4),xdata)-ydata);
sol = fminsearch(obj_fun, [0.5,1e-4,0.5,1e-4]);
a = sol(1);
b = sol(2)
c = sol(3);
d = sol(4)

xfit = linspace(min(xdata),max(xdata),100);
yfit = f(a, b, c, d, xfit);

figure(1);
semilogy(xdata, ydata, '+', 'MarkerSize', 10, 'LineWidth', 2)
hold on
semilogy(xfit, yfit, '-');

%% 3 parameters fminsearch optimization
d = 1e-3;
f = @(a,b,c,x) a.*exp(-b.*x) + c.*exp(-d.*x);
obj_fun = @(params) norm(f(params(1), params(2), params(3),xdata)-ydata);
sol = fminsearch(obj_fun, [0.5,1e-4,0.5,1e-4]);
a = sol(1);
b = sol(2)
c = sol(3);

xfit = linspace(min(xdata),max(xdata),100);
yfit = f(a, b, c, xfit);

figure(2);
semilogy(xdata, ydata, '+', 'MarkerSize', 10, 'LineWidth', 2)
hold on
semilogy(xfit, yfit, '-');

[Mathieu NOE]

Update

the problem can be even reduced to find 2 parameters , a and b for this equation :

d = 1e-3 (forced).

y = a.*exp(-b.*x) + (1-a).*exp(-d.*x);

this force the first point to have y = 1 ( y = a + 1 -a)

see the third plot :

clear; clc; clf; close all;
D_0 = [0.5 0.00001 0.5 0.00001]; %initial guess
xdata = [9.85	32.05	66.83	114.44	174.55	247.57	333.00	431.44	542.19	666.04	802.12	951.39	1113.38	1287.47	1474.87	1674.29	1887.11	2600.14	2863.78	3139.17	3428.23	4043.41	4369.45	4709.33	5425.99	5802.67	6193.38	6595.39];
ydata = [1	0.9569	0.9528	0.8894	0.8387	0.8995	0.7911	0.773	0.7523	0.7155	0.7086	0.6478	0.6269	0.6175	0.574	0.551	0.4991	0.4559	0.4449	0.4314	0.4212	0.407	0.3856	0.3511	0.3526	0.303	0.3148	0.2912];

%% 4 parameters fminsearch optimization
f = @(a,b,c,d,x) a.*exp(-b.*x) + c.*exp(-d.*x);
obj_fun = @(params) norm(f(params(1), params(2), params(3), params(4),xdata)-ydata);
sol = fminsearch(obj_fun, [0.5,1e-4,0.5,1e-4]);
a = sol(1);
b = sol(2)
c = sol(3);
d = sol(4)

xfit = linspace(min(xdata),max(xdata),100);
yfit = f(a, b, c, d, xfit);

figure(1);
semilogy(xdata, ydata, '+', 'MarkerSize', 10, 'LineWidth', 2)
hold on
semilogy(xfit, yfit, '-');

%% 3 parameters fminsearch optimization
d = 1e-3;
f = @(a,b,c,x) a.*exp(-b.*x) + c.*exp(-d.*x);
obj_fun = @(params) norm(f(params(1), params(2), params(3),xdata)-ydata);
sol = fminsearch(obj_fun, [0.5,1e-4,0.5]);
a = sol(1);
b = sol(2)
c = sol(3);

xfit = linspace(min(xdata),max(xdata),100);
yfit = f(a, b, c, xfit);

figure(2);
semilogy(xdata, ydata, '+', 'MarkerSize', 10, 'LineWidth', 2)
hold on
semilogy(xfit, yfit, '-');

%% 2 parameters fminsearch optimization
d = 1e-3;
f = @(a,b,x) a.*exp(-b.*x) + (1-a).*exp(-d.*x);
obj_fun = @(params) norm(f(params(1), params(2),xdata)-ydata);
sol = fminsearch(obj_fun, [0.5,1e-4]);
a = sol(1);
b = sol(2)

xfit = linspace(min(xdata),max(xdata),100);
yfit = f(a, b, xfit);

figure(3);
semilogy(xdata, ydata, '+', 'MarkerSize', 10, 'LineWidth', 2)
hold on
semilogy(xfit, yfit, '-');

[Alfredo Scigliani]

Forcing it to be 0.001 is still a problem for me, the answer should be around 0.0001 or even smaller. What I am looking for is to constraint it to be below it. Forcing it to a fixed value means assuming it to be that value.

[Mathieu NOE]

ok

I tried to see if my R² parameters evolves in a certain amount if I do a for loop to test with diffrent d values. here i test with 1000 values of d in log spacing from 10^-5 to 10^-2

there are two optimal points but in fact it's more or less the same solution , becuase parameters can be flipped (like b and d)

a = 3.914218228252134e-01 6.020081918678649e-01

b = 1.479928472725444e-03 1.073153563776005e-04

d = 1.079028791516184e-04 1.403289084785873e-03

Rsquared = 9.935951065959850e-01 9.938844449832044e-01

clear; clc; clf; close all;
xdata = [9.85	32.05	66.83	114.44	174.55	247.57	333.00	431.44	542.19	666.04	802.12	951.39	1113.38	1287.47	1474.87	1674.29	1887.11	2600.14	2863.78	3139.17	3428.23	4043.41	4369.45	4709.33	5425.99	5802.67	6193.38	6595.39];
ydata = [1	0.9569	0.9528	0.8894	0.8387	0.8995	0.7911	0.773	0.7523	0.7155	0.7086	0.6478	0.6269	0.6175	0.574	0.551	0.4991	0.4559	0.4449	0.4314	0.4212	0.407	0.3856	0.3511	0.3526	0.303	0.3148	0.2912];

%% 2 parameters fminsearch optimization
dd = logspace(-5,-2,1000);
for ci =1:numel(dd)
    d = dd(ci);
    f = @(a,b,x) a.*exp(-b.*x) + (1-a).*exp(-d.*x);
    obj_fun = @(params) norm(f(params(1), params(2),xdata)-ydata);
    sol = fminsearch(obj_fun, [0.5,1e-3]);
    a(ci) = sol(1);
    b(ci)  = sol(2);

    xfit = linspace(min(xdata),max(xdata),100);
    yfit = f(a(ci), b(ci), xfit);

    Rsquared(ci) = my_Rsquared_coeff(interp1(xdata,ydata,xfit),yfit); % correlation coefficient
end

% finding the best (or two best sets)
Rsquared_s = smoothdata(Rsquared,'gaussian',25);
Rsquared_s(Rsquared_s<=0.9) = 0.9;
ind = find(islocalmax(Rsquared_s));
figure(1);
semilogx(dd,Rsquared,dd(ind),Rsquared(ind),'dr');
title('R² vs d parameter ');
xlabel(' d parameter ');
ylabel(' R² ');

format long e
    a = a(ind)
    b  = b(ind)
    d  = dd(ind)
    Rsquared = Rsquared(ind)

    yfit1 = a(1).*exp(-b(1).*xfit) + (1-a(1)).*exp(-d(1).*xfit);
    yfit2 = a(2).*exp(-b(2).*xfit) + (1-a(2)).*exp(-d(2).*xfit);
    
figure(3);
semilogy(xdata, ydata, '+', 'MarkerSize', 10, 'LineWidth', 2)
hold on
semilogy(xfit, yfit1, '-', xfit, yfit2, '-');
title('Data fit ');
legend('data',['fit #1 , R² = ' num2str(Rsquared(1))],['fit #2 , R² = ' num2str(Rsquared(2))]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Rsquared = my_Rsquared_coeff(data,data_fit)
    % R2 correlation coefficient computation
    
    % The total sum of squares
    sum_of_squares = sum((data-mean(data)).^2);
    
    % The sum of squares of residuals, also called the residual sum of squares:
    sum_of_squares_of_residuals = sum((data-data_fit).^2);
    
    % definition of the coefficient of correlation is
    Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
    
    
end