Optimizing experimental data to theoretical solution-two parameters with infinite series

Hello, I am working on a problem that involves two constants in a non-linear equation, for which I would like to optimize the constants’ values based on experimental data. The equation used for optimization is

E_calc = 1 - 2 * (x(2)*a/x(1))^2 * sum [exp(-root_n^2 * x(1)*t/a^2)/(root_n^2 * (root_n^2 + (x(2)*a/x(1))^2 + x(2)*a/x(1))]

Where I would like to find the optimal values of x(1) and x(2). The objective function I would like to minimize is: (E_exp - E_calc)^2/N Where E_exp is the experimentally-measured value at time t, and N is the total number of points in time. I would like to implement some constraints on x(1) and x(2) and found the function fminsearchbnd on the file exchange, which performs the same operations as fminsearch while allowing to set bounds on the variables of interest. I am implementing this by providing educated initial guesses for x(1) and x(2), calculating the roots, and supplying the necessary values to a function (objectivefun). Two iterations occur so that the roots may be updated based on the optimal solution of x(1) and x(2) in the first iteration to arrive at a final answer in the second iteration. However, I have not been able to get the optimization procedure to work correctly. I have created data that I know the answer for and set the initial guesses to be very close, but no luck. I have also adjusted the tolerance with no luck. I suspect that I may be executing the objectivefun incorrectly. Does anyone have any guidance? The script and function that I have been using are shown below. I have used this example for guidance https://www.mathworks.com/help/matlab/examples/optimal-fit-of-a-nonlinear-function.html?s_tid=gn_loc_drop Thank you very much in advance!

%Known
x(1) = 4.87 * 10^-11;                           %Initial guess for D [m^2/s]
x(2) = 2.0 * 10^-7;                           %Initial guess for S [m/s]
a = 0.00159;

%Read in 'experimental' data
data = xlsread('input data');
t = data(:,1);
E_exp = data(:,2);
%bounds for x(1) and x(2)
lb = [1*10^-11,1.99*10^-7];  ub = [5.99*10^-11,2.01*10^-7];

for m = 1:2
x0 = [x(1);x(2)];

%Calculate roots
b = 1;
N = 10000;                                  %number of intervals
beta_max = 30;                              %Maximum beta
j = 0:N;                                    %'x' vector
beta_j = j/N*beta_max;                      %beta(x)
fun = (x(2)*a)/x(1) - beta_j.*tan(beta_j);              %Characteristic function over range of beta values
obj_fun = @(beta_j) ((x(2)*a)/x(1) - beta_j.*tan(beta_j));
for i = 1:N-1
    if fun(i) * fun(i+1) <= 0
        r_b(b,1) = b;
        r_b(b,2) = fzero(obj_fun,beta_j(i));
        b = b+1;
    else
    b = b;
    end
end
roots = zeros(length(r_b),2);
roots(:,1) = r_b(:,2);
%Run optimization procedure
[x, fval] = fminsearchbnd((@(x) objectivefcn(x,t,E_exp,roots(:,1),a)),x0,lb,ub);

End

function f = objectivefcn(x,t,E_exp,roots,a)
%Computing objective function to optimize D/s
%   x(1) = D [m^2/s]
%   x(2) = S [m/s]
%   t = time [s]
%   E_exp = experimental data
%   beta_n
%   a = length in transport direction [m]
A = zeros(length(t),1);

   for k = 1:length(t)
       for z = 1:length(roots)
       E_T1(z,k) = exp(-roots(z,1)^2 * (x(1)*t(k))/a^2) / ((roots(z,1)^2) * (roots(z,1)^2 + (x(2)*a/x(1))^2 + (x(2)*a/x(1))));
       end
       A(k,1) = 1 - (2 * ((x(2)*a/x(1)))* sum(E_T1(:,k)));

   end
   C = A\E_exp;
   Z = A*C;
   f = norm(Z-E_exp);
   end

[Torsten]

1. Since x(1) and x(2) change during the optimization and if the roots are defined as you calculate them, they will also change during the optimization. So you can't calculate the roots in advance. But I doubt that the roots are defined the way you calculate them.

2. I don't understand the last three lines in "objectivefcn". If C=A\E_Exp, then A*C=E_exp, thus Z=E_exp and f=norm(Z-E_exp)=0.

Best wishes

Torsten.

[Daniel Way]

Torsten, Thank you for pointing that out.

Answer 1

[Matt J]

The objective function I would like to minimize is: (E_exp - E_calc)^2/N

If so, then shouldn't your objective be this,

function f = objectivefcn(x,t,E_exp,roots,a)

   N=length(t);

    A = zeros(N,1);

    rsq=roots(:,1).^2;

    denom= (rsq .* ( rsq + (  x(2)*a/x(1) )^2 + (x(2)*a/x(1))  )  );

     for k = 1:N

         E_T1 = exp(-rsq * (x(1)*t(k))/a^2) ./denom;

         A(k) = 1 - (2 * ((x(2)*a/x(1)))* sum(E_T1);

     end

     f = norm(A(:)-E_exp(:)).^2/N;

end

[Daniel Way]

That did the trick. Thank you very much.

[Torsten]

@Daniel Way:

According to your image file,

A(k) = 1 - 2 * (x(2)*a/x(1))^2* sum(E_T1);

instead of

A(k) = 1 - (2 * ((x(2)*a/x(1)))* sum(E_T1);

Best wishes

Torsten.