24 views (last 30 days)

Show older comments

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!

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

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.

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

Start Hunting!