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

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)

[Ameer Hamza]

Can you attach your dataset so that it will be easy for us to suggest a solution?

[Arthur Leonard]

% my code is the the following and the dataset is attached, ty for your time!
function devoir()
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
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)
end

Answer 1

[Ameer Hamza]

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

[Arthur Leonard]

In fact, I don't have access to ga() function but your second solution works perfectly thank you very much for your help and your time!

[Ameer Hamza]

Glad to be of help

Answer 2

[Image Analyst]

My guess would be to just optimize the parameter as a floating point value first. Then take the integers on either side of that value and compute the residuals to figure out which integer gives the closest overall to the data.

[Arthur Leonard]

Thank you for your quick and helpful answer. I guess it's a good idea as the first step should be easy to do with this dataset. I just have one more question, is there a Matlab function I can use to compute the residuals?