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

fminsearchfminsearchbndparameter estimation/optimization

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

Best Answer

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
Related Question