MATLAB: Hot to fit two curves under interdependent constraint

Question

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

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