MATLAB Answers

Hot to fit two curves under interdependent constraint ?

24 views (last 30 days)
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!

Accepted Answer

Alan Stevens
Alan Stevens on 22 Feb 2021
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
  2 Comments
Alan Stevens
Alan Stevens on 26 Feb 2021
"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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!