Multivariable nonlinear curve fitting with equality and inequality constraints

Dear MatLab experts,

currently I'm facing small problems implementing a curve fitting and hoping to find your assistance.

Given the vectors u and v containing the measured data points (u(i), v(i)), I want to fit the polynominal function v = a*u^ea + b*u^eb + c*u^ec as good as possible. a, b, c, as well as ea, eb and ec are the degrees of freedom. Unfortunately there are the following constraints, which drive me crazy:

• a, b, c, ea, eb and ec have to be positive
• a, b and c have to sum up to 1

During my research I developed a sketch to find a fit for one data pair:

Objective function definition:

function fval = objectiveFunction(x)
  % Fit to the following mathematical equation (example with u = 1, v = 4)
  % v = a*u^ea + b*u^eb + c*u^ec
  fval = x(1)*1^x(4) + x(2)*1^x(5) + x(3)*1^x(6) - 4;
end

My attempt to find the optimum parameters of the function:

%%--- Define constraints ---
% There is no inequality constraint
A = [];
b= [];

% Equality constraint: a, b and c have to sum up to 1
Aeq = [ 1 1 1 0 0 0 ];
beq = 1;

%%--- Set upper and lower bounds ---
lb = zeros(1,6);                    % a b c ea eb ec
ub = [ Inf Inf Inf Inf Inf Inf ];   % a b c ea eb ec

%%--- Initial guess ---
x0 = zeros(6,1);

%%--- Find the optimum ---
x = fmincon( @(x)objectiveFunction(x), x0, A, b, Aeq, beq, lb, ub );

But I don't have any clue how to efficiently extend this code for a arbitrary number of data points. My first thought is to adjust the objective function to

function fval = objectiveFunction(x)
  % Minimize deviation for each measured value
  % (a*u(1)^ea + b*u(1)^eb + c*u(1)^ec - v(1))^2 + (a*u(2)^ea + b*u(2)^eb + c*u(2)^ec - v(2))^2 + ...
  fval = (x(1)*x(7)^x(4) + x(2)*x(7)^x(5) + x(3)*x(7)^x(6) - x(8))^2; % + ...
end

and subscribe the new degrees of freedom to the equality constraints, so that x(7) = u(1), x(8) = v(1) and so on.

From my point of view it shouldn't be that complicated to extend the matrices Aeq and beq inside of a while-loop. A fortiori I have problems figuring out how to extend the objective function efficiently.

I would be very greatful, if you could give me a hint on either my approach or on how to get my curve fit in a more advantageous way. Many thanks in advance!

Best regards, Chris

Answer 1

[Torsten]

Usually in curve fitting, the objective function is given by

fval = sum((x(4)*u.^x(1)+x(5)*u.^x(2)+x(6)*u.^x(3)-v).^2)

Best wishes

Torsten.

[Chris Kittl]

Thanks for your solution, Torsten! You gave me the right idea to solve my problem!

For anyone googling this and wondering what could lead the way:

Test data:

u = linspace(1, 3, 3);
v = u.^2;

Code for finding the optimized fit:

% This script fits a curve to a set of data pairs while complying different
% constraints.
% The objective function is v = a*u^ea + b*u^eb + c*u^ec
% Constraint 1: a, b, c, ea, eb and ec have to be >= 0
% Constraint 2: a + b + c = 1

%%--- Defining constraints, boundaries and other stuff ---
% There is no inequality constraint
A = [];
b= [];

% Equality constraint: a, b and c have to sum up to 1
% Other euqality constraints are set up by the sample points
Aeq = [1 1 1 0 0 0];
beq = 1;

% Set upper and lower bounds
lb = zeros(1,6);    % a b c ea eb ec
ub = Inf(1,6);      % a b c ea eb ec

% Initial guess
x0 = rand(6,1);

%%--- Find the optimum ---
options = optimoptions( @fmincon, 'TolFun', 1E-12 );
x = fmincon( @(x)sum((x(1)*u.^x(4)+x(2)*u.^x(5)+x(3)*u.^x(6)-v).^2), ...
  x0, A, b, Aeq, beq, lb, ub, [], options );
a = x(1);
b = x(2);
c = x(3);
ea = x(4);
eb = x(5);
ec = x(6);

%%--- Check the goodness ---
vtest = a*u.^ea + b*u.^eb + c*u.^ec;
RMSE = sqrt(mean((vtest - v).^2));

plot( u, v, '*');
hold on;
plot( u, vtest);
hold off;
title(sprintf('Goodness check (RMSE = %f)',RMSE))
xlabel('x');
ylabel('y');

Most probably this isn't the most efficient solution, achieving the best fit, but it works fine for my objective. Although the problem is solved, please feel free to optimize my approach.

Many thanks for your quick assistance, Torsten!

Best regards, Chris