MATLAB: How to perform a non linear curve fitting with the parameter to optimize being an integer

Question

I need to perform a curve fitting on the parameter n in the following code but don't know how to do as my parameter n is an integer.

I tried to simply perform the curve fitting and round the result to the nearest decimal but it's not possible as their is a factorial in my function. Could anyone help me with this one? I also tried this:

fact = factorial (round(n_reac)-1);

Any help would be welcome

Here is my Matlab code:

function E_n = fun (n_reac,time)

tau = 1.2584

tau_series = tau/n_reac;

fact = factorial (n_reac-1);

E_n = (time.^(n_reac-1).*exp((-time)./tau_series))./(fact.*(tau_series^n_reac));

end

n_optimal = lsqcurvefit(@fun, 1, time, esp)

Answer 1

First, note that you are using the factorial function in your model. So you cannot simply use the factorial function. In Matlab it is only defined for non-negative inputs. Instead use gamma() function. It is equivalent to factorial(n) = gamma(n+1) for integer value, but it also extended to non-integer and negative values.

Now, for solving the problem, if you have the global optimization toolbox, then my first suggestion would be to use genetic algorithm ga() to solve the problem. ga() has the option to enforce integer constraints. Following is the code

Table = csvread('tracer_data_reactorA_new.csv')';
time = Table(1,:);
tracer_conc = Table(2,:);
int_tot = 0;
for i = 1:size(time,2)-1
  int_tot = int_tot + (0.5*tracer_conc(i)+0.5*tracer_conc(i+1))*(time(i+1)-time(i));
end
esp = tracer_conc/int_tot;
obj_fun = @(n_reac) norm(fun(n_reac, time) - esp);
n_reac_sol = ga(obj_fun, 1, [], [], [], [], [], [], [], 1);
function E_n = fun(n_reac,time)
    tau = 1.2584;
    tau_series = tau/n_reac;
    fact = gamma(n_reac);
    E_n = (time.^(n_reac-1).*exp((-time)./tau_series))./(fact.*(tau_series^n_reac));
end

If you don't have ga() function available to you, then you can follow the Image Analyst's advice and solve the problem for real numbers and finally round off the value. Note that this may not always work, since the objective function may not always be convex. However, In this case, it seems to work and both ga and lsqcurvefit gives very close output

Table = csvread('tracer_data_reactorA_new.csv')';
time = Table(1,:);
tracer_conc = Table(2,:);
int_tot = 0;
for i = 1:size(time,2)-1
  int_tot = int_tot + (0.5*tracer_conc(i)+0.5*tracer_conc(i+1))*(time(i+1)-time(i));
end
esp = tracer_conc/int_tot;
obj_fun = @(n_reac) norm(fun(n_reac, time) - esp);
n_optimal = lsqcurvefit(@fun, 1, time, esp);
n_optimal = round(n_optimal); % you may also check the value of objective function
                              %  at round(n_optimal) and round(n_optimal)+1 round(n_optimal)-1
function E_n = fun(n_reac,time)
    tau = 1.2584;
    tau_series = tau/n_reac;
    fact = gamma(n_reac);
    E_n = (time.^(n_reac-1).*exp((-time)./tau_series))./(fact.*(tau_series^n_reac));
end