Lsqcurvefit, upper bounds that depend on other upper bounds

Question

In trying to keep the question simple, I'm using MATLAB and trying to create an upper bound for lsqcurvefit to use where the sum of my upper boundaries on the parameters cannot exceed 100.

Something like

lb = [0, 0, 0]; ub = [a, b, c];

where a+b+c < 100

Is there a simple elegant way to pass that type of information into the lsqcurvefit command?

I've tried setting limiters in the ode45 portion with if statements, I've tried 100s in each of the upper bound spots to see if the lsq came out naturally following that limitation. What I'm expecting is that the parameters I'm getting out of lsqcurvefit should sum to something less than or equal to 100.

Answer 1

The lsqcurvefit function doesn't have inputs for constraints between optimisation parameters as far as I can tell, so you'll likely need to use a more generic optimisation function like fmincon.

Here is a comparison, where I take a quadratic function, which accepts a 3-element array of coefficients abc and axis points x:

fun = @(abc, x) abc(1)*x.^2 + abc(2)*x + abc(3);

We can evaluate the true values at these points:

abc = [2, 5, 3]; % note the true sum of these coeffs is 10!
y = fun(abc, x);

And set up some constraints for out optimisation:

abcMax = 9;     % The max permitted sum of the coefficients abc
abc0 = [1 2 3]; % Initial guess for abc during optimisation

ub = [abcMax,abcMax,abcMax]; % Uppper bounds for abc, use the max sum of each
lb = [0,0,0];                % Lower bounds for abc, using 0 for each

Using lsqcurvefit, you can only ensure that each element of abc is within the bounds, using something like

abcLSQ = lsqcurvefit( fun, abc0, x, y, lb, ub );

However, we can construct an additional inequality constraint for fmincon and use that instead:

% Minimise the absolute error F(abc) = ( fun(abc,x) - y ).^2
% For scalar minimisation, we minimise sum(F)
% such that A*abc.' <= B
% and lb <= abc <= ub
A = [1 1 1];  % 1a + 1b + 1c
B = abcMax; % a + b + c <= abcTotal
abcFMC = fmincon( @(abc) sum((fun(abc,x) - y).^2), abc0.', A, B, [], [], lb, ub );

We can do some plotting to confirm the result - I've added the output coefficient arrays to the legend

% Create an x vector for fitting which is slightly longer than x
xfit = linspace(x(1)-1,x(end)+1,2*numel(x));
% Fit using the two coefficient sets
yfitLSQ = fun( abcLSQ, xfit );
yfitFMC = fun( abcFMC, xfit );

% Plot
figure(1); clf; 
hold on;
fstr = @(str, abc) sprintf( '%s, [%.2f, %.2f, %.2f]', str, abc(1), abc(2), abc(3) );
plot( x, y, 'ob', 'displayname', fstr('input data', abc), 'markersize', 5 );
plot( xfit, yfitLSQ, 'r', 'displayname', fstr('lsqcurvefit', abcLSQ), 'linewidth', 1 );
plot( xfit, yfitFMC, 'k--', 'displayname', fstr('fmincon', abcFMC), 'linewidth', 1 );
legend( 'location', 'northwest' ); grid on;
title( 'Fits for equation of form a*x^2 + b*x + c' );