Hot to fit two curves under interdependent constraint ?
1 view (last 30 days)
Show older comments
Stevan Pecic
on 22 Feb 2021
Commented: Alan Stevens
on 26 Feb 2021
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:
- 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.
- 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!
0 Comments
Accepted Answer
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
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.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!