Fit theoretical x-y data set to experimental data set

I have a function to model the force-displacement relation between an interface and a probe. The function outputs two data sets (the probe movement and corresponding force) based on three inputs. I have experimental data and essentially want to find what value for the inputs results in the closest fit according to my model. The problem looks like this,

function [f,d] = myFunc(x)

where, f and d are arrays of data points that I want to fit as closely as possible to the experimental data i have, [f_exp, d_exp].

The crux of the problem to me is I am not sure how to constrain my outputs in order to data mine the inputs. I have so far been trying to use the genetic algorithm function in the optimisation tool box but am not sure if this is the correct approach.

Thanks for any help.

Cheers,

Will

[Star Strider]

‘... constrain my outputs in order to data mine the inputs.’

Please go into a bit more detail on this. Do you just want to do an unconstrained optimisation to get a good fit to your data, or are your parameters (that I assume are ‘x’) themselves constrained in some way?

[Will Meg]

As a starting point I would just like unconstrained optimisation to get a good fit to my data. My parameter 'x' consists of three parameters x(1), x(2) and x(3), which are not constrained. But want the output data sets to be as close a match as possible to the experimental data i have obtained.

[Image Analyst]

Will, where are the screenshots? Where did you attach the data? Why not make it easy for people to answer you, not hard?

If you want it to fit your experimental data perfectly, then use spline interpolation. If you want a regression, then there are several ways but you might have to pick a model equation (which we have no guesses at since you didn't attach a screenshot).

[Will Meg]

So my function solves the augmented Young-Laplace equation using a 4th order runga kutta numerical method, and based on this solution I can get force displacement data, as shown in my script below;

function [FX0] = RK4AL_func(myin)
%% Augmented Young Laplace
counter = 0;
for d0 = 160:1.5:1558
hold off
counter = counter+1;
A=4e-20;
e=80;                                               % Relative permitivity water
e0=8.85e-12;                                        % Permitivity of free space
k=100;                                               % debeye length nm
dRho = 150;                                         % delta density kg/m^3 (same as mg/cm^3)
zp1=myin(1);                                          % zeta potential 1
zp2=myin(2);                                          % zeta potential 2
%d0=195;                                            % minimum interfacial gap nm
a=10e-6;                                            % bead radius nm
g=myin(3);                                            % surface tension N/m
h=0.1;                                              % step size
x = 0.1:h:25;                                       % Calculates upto y(1)
y = zeros(1,length(x)); 
z = zeros(1,length(x)); 
y(1) = d0;                                          % initial condition
z(1) = -0.000000000000001;                          % initial condition
% F_xy = @(t,r) 3.*exp(-t)-0.4*r;                   % change the function as you desire
F_xyz = @(x,y,z) z;                                 % change the function as you desire
%G_xyz = @(x,y,z) (2-((20*(((80*8.85e-12*(2*-0.02*-0.02*(exp(-1*y/96.2)+exp(-3*y/96.2))+(((-0.02)^2+(-0.02)^2)*exp(-2*y/96.2))))/((96.2^2)*((1-exp(-2*y/96.2))^2)))*1e18)-(4e-20*1e27/(6*pi*(y^3))))/50))-(z/x);  
G_xyz = @(x,y,z) d0*(2-((a*(((e*e0*(2*zp1*zp2*(exp(-1*y/k)+exp(-3*y/k))+(((zp1^2)+(zp2^2))*exp(-2*y/k))))/((k^2)*((1-exp(-2*y/k))^2)))*1e18)-(A*1e27/(6*pi*(y^3))))/g))-(z/x);
%DP(i) = (((e*e0*(2*zp1*zp2*(exp(-1*y(i)/k)+exp(-3*y(i)/k))+((zp1^2+zp2^2)*exp(-2*y(i)/k))))/((k^2)*((1-exp(-2*y(i)/k))^2)))*1e18)-(A*1e27/(6*pi*(y(i)^3)));
for i=1:(length(x)-1)                               % calculation loop
    k_1 = F_xyz(x(i),y(i),z(i));
    L_1 = G_xyz(x(i),y(i),z(i));
    k_2 = F_xyz(x(i)+0.5*h,y(i)+0.5*h*k_1,z(i)+0.5*h*L_1);
    L_2 = G_xyz(x(i)+0.5*h,y(i)+0.5*h*k_1,z(i)+0.5*h*L_1);
    k_3 = F_xyz((x(i)+0.5*h),(y(i)+0.5*h*k_2),(z(i)+0.5*h*L_2));
    L_3 = G_xyz((x(i)+0.5*h),(y(i)+0.5*h*k_2),(z(i)+0.5*h*L_2));
    k_4 = F_xyz((x(i)+h),(y(i)+k_3*h),(z(i)+L_3*h)); % Corrected        
    L_4 = G_xyz((x(i)+h),(y(i)+k_3*h),(z(i)+L_3*h));

    y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  % main equation
    z(i+1) = z(i) + (1/6)*(L_1+2*L_2+2*L_3+L_4)*h;  % main equation

end

%% Drop Profile

hold on
DrPr = (y-d0)/d0;                                   %DrPr is the drop profile (D(t)-D0)/D0
plot(x,DrPr)

%% Finding G, H, F and X 

for j = 1:length(x)
DP(j) = (((e*e0*(2*zp1*zp2*(exp(-1*y(j)/k)+exp(-3*y(j)/k))+((zp1^2+zp2^2)*exp(-2*y(j)/k))))/((k^2)*((1-exp(-2*y(j)/k))^2)))*1e18)-(A*1e27/(6*pi*(y(j)^3)));
dG(j) = a*d0*x(j)*DP(j)/g;
dH(j) = a*d0*x(j)*log(x(j))*DP(j)/g;
end

for jj = 2:length(x);
G(jj) = (dG(jj)+dG(jj-1))*h/2;
H(jj) = (dH(jj)+dH(jj-1))*h/2;
end

Gsum(counter) = sum(G);
Hsum(counter) = sum(H);

B = 0.57721566 + log(sqrt(a)/2*sqrt(g/(dRho*9.81)));
FX0(counter,1) = (d0 - Hsum(counter) - Gsum(counter)*((1/2*log(d0))+B))-160; %Piezo Movement, X0
FX0(counter,2) = (2*pi*g*Gsum(counter)/1000)*1e6; %Force observed by AFM, F
end

end

I then have experimental force displacement data that I would like to get as close a fit as possible to by varying three inputs to my function. The force - displacement data is essentially just x-y data and looks like this.

So I want to am trying to use the optimization function so that the exported data from my funciton matches this data as closely as possible.

Answer 1

[Star Strider]

I have no idea what your actual function is, or what you are doing.

I would do something like this:

function fd = myFunc(x)
fd(:,1) = ...;                                                  % Approximate ‘f_exp’
fd(:,2) = ...;                                                  % Approximate ‘d_exp’
end

fitnessfcn = @(x) norm([f_exp(:), d_exp(:)] - myFunc(x));

Ideally, if you are fitting data, ‘myFunc’ should be a function of your independent variable and the parameter vector. I have no idea how ‘f_exp’ relates to ‘d_exp’, or to your parameters, so this is illustrative only.

EDIT —

Since you’re fitting differential equations to data, see:

Monod kinetics and curve fitting (link)

Solving Coupled Differential Equations (link)

Parameter Estimation for a System of Differential Equations (link)