Hot to fit two curves under interdependent constraint ?

Hello,

I'm trying to fit two functions ; first one is the calibration curve y(z) and the second one is the profile of the measured quantity z(x). This calibration curve is modeled as polynomial:

and profile is modeled as linear function:

Measurement results to which these two functions must be fitted are three vectors y1(x),y2(x),y3(x) , which are related to each other as:

My questions:

1. How to impose these relations as a constraint to a least-squares/optimization problem ? So far I belive that fmincon is the appropriate tool for the job, but i'm not sure of the exact implementation.
2. Should i loop over the given measurment vectors and averege out the parameters or is there a method of finding unique parameters a,b,k and m for these vectors by default ?

One possible approach:

Given that these vectors depened on x, it seems natural to combine the two functions and obtain the y(x) functional form. This way, problem is reduced to fitting the function:

I've tried implementing this as writing the objective function as:

separating k and m varables for each measurement vector and imposing the constraint to these variables as:

;

;

Thereby, the goal is to find a,b,k1 and m1, as these correspond to the parameters k and m.

This has been implemented as following:

x = linspace(1,100,100);
% Measurement results:
y1=[0.07469056,0.07378624,0.07288704,0.07199296,0.071104,0.07022016,0.06934144,0.06846784,0.06759936,0.066736,0.06587776,0.06502464,0.06417664,0.06333376,0.062496,0.06166336,0.06083584,0.06001344,0.05919616,0.058384,0.05757696,0.05677504,0.05597824,0.05518656,0.0544,0.05361856,0.05284224,0.052071040,0.051304960,0.05054400,0.04978816,0.04903744,0.04829184,0.04755136,0.046816,0.04608576,0.04536064,0.04464064,0.04392576,0.043216,0.04251136,0.04181184,0.04111744,0.04042816,0.039744,0.03906496,0.03839104,0.03772224,0.03705856,0.0364,0.03574656,0.03509824,0.03445504,0.03381696,0.033184,0.03255616,0.03193344,0.03131584,0.03070336,0.030096,0.02949376,0.02889664,0.02830464,0.02771776,0.027136,0.02655936,0.02598784,0.02542144,0.02486016,0.024304,0.02375296,0.02320704,0.02266624,0.02213056,0.0216,0.02107456,0.02055424,0.02003904,0.01952896,0.019024,0.01852416,0.01802944,0.01753984,0.01705536,0.016576,0.01610176,0.01563264,0.01516864,0.01470976,0.014256,0.01380736,0.01336384,0.01292544,0.01249216,0.012064,0.01164096,0.01122304,0.01081024,0.01040256,0.01];
y2=[0.23624224,0.23310496,0.22998816,0.22689184,0.223816,0.22076064,0.21772576,0.21471136,0.21171744,0.208744,0.20579104,0.20285856,0.19994656,0.19705504,0.194184,0.19133344,0.18850336,0.18569376,0.18290464,0.180136,0.17738784,0.17466016,0.17195296,0.16926624,0.1666,0.16395424,0.16132896,0.15872416,0.15613984,0.153576,0.15103264,0.14850976,0.14600736,0.14352544,0.141064,0.13862304,0.13620256,0.13380256,0.13142304,0.129064,0.12672544,0.12440736,0.12210976,0.11983264,0.117576,0.11533984,0.11312416,0.11092896,0.10875424,0.1066,0.10446624,0.10235296,0.10026016,0.09818784,0.096136,0.09410464,0.092093760,0.09010336,0.0881334400,0.0861840,0.084255040,0.08234656,0.080458560,0.078591040,0.0767440,0.07491744,0.07311136,0.07132576,0.06956064,0.067816,0.06609184,0.06438816,0.06270496,0.06104224,0.0594,0.05777824,0.05617696,0.05459616,0.05303584,0.051496,0.049976640,0.04847776,0.04699936,0.04554144,0.044104,0.04268704,0.04129056,0.03991456,0.03855904,0.037224,0.03590944,0.03461536,0.03334176,0.03208864,0.030856,0.02964384,0.02845216,0.02728096,0.02613024,0.025];
y3=[0.48465504,0.477956160,0.47130336,0.46469664,0.458136,0.45162144,0.44515296,0.43873056,0.43235424,0.426024,0.41973984,0.41350176,0.40730976,0.40116384,0.395064,0.38901024,0.38300256,0.37704096,0.37112544,0.365256,0.35943264,0.35365536,0.34792416,0.34223904,0.3366,0.33100704,0.32546016,0.31995936,0.31450464,0.309096,0.30373344,0.29841696,0.29314656,0.28792224,0.282744,0.27761184,0.27252576,0.26748576,0.26249184,0.257544,0.25264224,0.24778656,0.24297696,0.23821344,0.233496,0.22882464,0.22419936,0.219620160,0.21508704,0.2106,0.20615904,0.20176416,0.19741536,0.193112640000000,0.188856,0.18464544,0.18048096,0.17636256,0.17229024,0.168264,0.16428384,0.16034976,0.15646176,0.15261984,0.148824,0.14507424,0.14137056,0.13771296,0.13410144,0.130536,0.12701664,0.12354336,0.12011616,0.11673504,0.1134,0.11011104,0.10686816,0.10367136,0.10052064,0.097416,0.09435744,0.09134496,0.08837856,0.08545824,0.082584,0.07975584,0.07697376,0.07423776,0.07154784,0.068904,0.06630624,0.06375456,0.06124896,0.05878944,0.056376,0.05400864,0.05168736,0.04941216,0.04718304,0.045];

% Point for minimization
xvalue = 50;
y1value = y1(xvalue);
y2value = y2(xvalue);
y3value = y3(xvalue);

% Objective function:
problem.objective = @(x)(y1value-(x(1)*(x(3))^2)*xvalue^2-(2*x(1)*x(3)*x(6)+x(2)*x(3))*xvalue-(x(2)*x(6)+x(1)*(x(6))^2))^2 + (y2value-(x(1)*(x(4))^2)*xvalue^2-(2*x(1)*x(4)*x(7)+x(2)*x(4))*xvalue-(x(2)*x(7)+x(1)*(x(7))^2))^2 + (y3value-(x(1)*(x(5))^2)*xvalue^2-(2*x(1)*x(5)*x(8)+x(2)*x(5))*xvalue-(x(2)*x(8)+x(1)*(x(8))^2))^2;
% x(1) = a, x(2) = b, x(3) = k1, x(4) = k2, x(5) = k3, x(6) = m1, x(7) = m2, x(8) = m3

%% Real parameter values: 
% x(1) = -4e-6
% x(2) = 3e-4
% x(3) = -0.8
% x(4) = -1.6
% x(5) = -2.4
% x(6) = 105
% x(7) = 210
% x(8) = 315
%%

problem.lb = [-1 -1 -5 -5 -5 0 0 0];   % Lower bound
problem.ub = [1 1 0 0 0 1000 1000 1000];  % Upper bound
problem.x0 = [-0.5,0.5,-2,-4,-6,100,200,300]; % Initial solution
problem.nonlcon = @film_constraint; % Constraints function handle

x = fmincon(problem)

function [c, ceq] = film_constraint(x)
    c = [];
    ceq = [x(7)-2*x(6); x(8)-3*x(6); x(5)-3*x(3);x(4)-2*x(3)];
end

So far it seems like this approach is not giving the plausible results.

I would be very thankful for any suggestion and/or comment!

Answer 1

[Alan Stevens]

How about just using fminsearch

x = linspace(1,100,100);
% Measurement results:
y1=[0.07469056,0.07378624,0.07288704,0.07199296,0.071104,0.07022016,0.06934144,0.06846784,0.06759936,0.066736,0.06587776,0.06502464,0.06417664,0.06333376,0.062496,0.06166336,0.06083584,0.06001344,0.05919616,0.058384,0.05757696,0.05677504,0.05597824,0.05518656,0.0544,0.05361856,0.05284224,0.052071040,0.051304960,0.05054400,0.04978816,0.04903744,0.04829184,0.04755136,0.046816,0.04608576,0.04536064,0.04464064,0.04392576,0.043216,0.04251136,0.04181184,0.04111744,0.04042816,0.039744,0.03906496,0.03839104,0.03772224,0.03705856,0.0364,0.03574656,0.03509824,0.03445504,0.03381696,0.033184,0.03255616,0.03193344,0.03131584,0.03070336,0.030096,0.02949376,0.02889664,0.02830464,0.02771776,0.027136,0.02655936,0.02598784,0.02542144,0.02486016,0.024304,0.02375296,0.02320704,0.02266624,0.02213056,0.0216,0.02107456,0.02055424,0.02003904,0.01952896,0.019024,0.01852416,0.01802944,0.01753984,0.01705536,0.016576,0.01610176,0.01563264,0.01516864,0.01470976,0.014256,0.01380736,0.01336384,0.01292544,0.01249216,0.012064,0.01164096,0.01122304,0.01081024,0.01040256,0.01];
y2=[0.23624224,0.23310496,0.22998816,0.22689184,0.223816,0.22076064,0.21772576,0.21471136,0.21171744,0.208744,0.20579104,0.20285856,0.19994656,0.19705504,0.194184,0.19133344,0.18850336,0.18569376,0.18290464,0.180136,0.17738784,0.17466016,0.17195296,0.16926624,0.1666,0.16395424,0.16132896,0.15872416,0.15613984,0.153576,0.15103264,0.14850976,0.14600736,0.14352544,0.141064,0.13862304,0.13620256,0.13380256,0.13142304,0.129064,0.12672544,0.12440736,0.12210976,0.11983264,0.117576,0.11533984,0.11312416,0.11092896,0.10875424,0.1066,0.10446624,0.10235296,0.10026016,0.09818784,0.096136,0.09410464,0.092093760,0.09010336,0.0881334400,0.0861840,0.084255040,0.08234656,0.080458560,0.078591040,0.0767440,0.07491744,0.07311136,0.07132576,0.06956064,0.067816,0.06609184,0.06438816,0.06270496,0.06104224,0.0594,0.05777824,0.05617696,0.05459616,0.05303584,0.051496,0.049976640,0.04847776,0.04699936,0.04554144,0.044104,0.04268704,0.04129056,0.03991456,0.03855904,0.037224,0.03590944,0.03461536,0.03334176,0.03208864,0.030856,0.02964384,0.02845216,0.02728096,0.02613024,0.025];
y3=[0.48465504,0.477956160,0.47130336,0.46469664,0.458136,0.45162144,0.44515296,0.43873056,0.43235424,0.426024,0.41973984,0.41350176,0.40730976,0.40116384,0.395064,0.38901024,0.38300256,0.37704096,0.37112544,0.365256,0.35943264,0.35365536,0.34792416,0.34223904,0.3366,0.33100704,0.32546016,0.31995936,0.31450464,0.309096,0.30373344,0.29841696,0.29314656,0.28792224,0.282744,0.27761184,0.27252576,0.26748576,0.26249184,0.257544,0.25264224,0.24778656,0.24297696,0.23821344,0.233496,0.22882464,0.22419936,0.219620160,0.21508704,0.2106,0.20615904,0.20176416,0.19741536,0.193112640000000,0.188856,0.18464544,0.18048096,0.17636256,0.17229024,0.168264,0.16428384,0.16034976,0.15646176,0.15261984,0.148824,0.14507424,0.14137056,0.13771296,0.13410144,0.130536,0.12701664,0.12354336,0.12011616,0.11673504,0.1134,0.11011104,0.10686816,0.10367136,0.10052064,0.097416,0.09435744,0.09134496,0.08837856,0.08545824,0.082584,0.07975584,0.07697376,0.07423776,0.07154784,0.068904,0.06630624,0.06375456,0.06124896,0.05878944,0.056376,0.05400864,0.05168736,0.04941216,0.04718304,0.045];

K0 = [0.1 0 0 1];   % initial guesses for [a b k m]
K = fminsearch(@(K) fn(K,y1,y2,y3,x), K0);

a = K(1); b = K(2); k = K(3); m = K(4);
disp([a b k m])
 Z1 = k*x + m;
 Y1 = a*Z1.^2 + b*Z1;
 Z2 = 2*Z1;
 Y2 = a*Z2.^2 + b*Z2;
 Z3 = 3*Z1;
 Y3 = a*Z3.^2 + b*Z3;

plot(x,y1,'ro',x,y2,'bo',x,y3,'ko'),grid
hold on
plot(x,Y1,'r',x,Y2,'b',x,Y3,'k')
xlabel('x'),ylabel('y')
legend('y1 data','y2 data','y3 data', 'y1 fit', 'y2 fit', 'y3 fit')
text(10,0.52,['[a b k m] = ',num2str([a b k m],4)]);

function F = fn(K,y1,y2,y3,x)

         a = K(1); b = K(2); k = K(3); m = K(4);
         
         Z1 = k*x + m;
         Y1 = a*Z1.^2 + b*Z1;
         Z2 = 2*Z1;
         Y2 = a*Z2.^2 + b*Z2;
         Z3 = 3*Z1;
         Y3 = a*Z3.^2 + b*Z3;
         
         F = norm(Y1-y1) + norm(Y2-y2) + norm(Y3-y3);
end

[Stevan Pecic]

Hi, thanks for the answer!

Your code does not seem to give results close to the exact values given. However, your formulation of the problem (definition of the objective function) seems to make much more sense than the previous one. After replacing fminsearch with fmincon (and trying different algorithms), i have arrived at something useful (results are somewhere within 10 percent of the exact values). Thank you for that improvement!

Now i'm stuck with the problem of forcing the algorithm not to stuck in the local minimum and find the global one. Any ideas on which optimization options to tweak in order to escape local minima ?

[Alan Stevens]

"Your code does not seem to give results close to the exact values given."

I'm surprised, as the fit looks excellent by eye!

To avoid getting stuck in a local minimum you could try a totally different sort of optimisation method, such as particle swarm optimisation (PSO). Search the File Exchange for useful routines.