How to optimise an objective function with a summation of integrals in Matlab?

Question

I'm trying to replicate the optimisation problem in the following paper:

Statistical modelling of directional wind speeds using mixtures of von Mises distributions: Case study

The objective function is as follows:

enter image description here

The code I have used to implement the whole process is shown below. I have the problem in the definition of the objective function

file = "dates.nc";
WD = ncread(file, 'WD');


WD_clean = WD(~isnan(WD));
WD_rad = deg2rad(WD_clean);

T= 360;

Limites_Sectores_T = deg2rad(0:360/T:360); % en radianes


Muestras_Sectores_T = histcounts(WD_rad, Limites_Sectores_T);


total_samples = numel(WD_rad);
P = cumsum(Muestras_Sectores_T) / total_samples; 


N = 4;

Limites_Sectores_N = deg2rad(0:(360/4):360);


Muestras_Sectores_N = histcounts(WD_rad, Limites_Sectores_N);

%s_j y c_j

%  s_j y c_j 
s_j = zeros(1, N);
c_j = zeros(1, N);
for j = 1:N
    
    s_j(j) = sum(sin(WD_rad(WD_rad >= Limites_Sectores_N(j) & WD_rad < Limites_Sectores_N(j+1)))) / Muestras_Sectores_N(j);
    c_j(j) = sum(cos(WD_rad(WD_rad >= Limites_Sectores_N(j) & WD_rad < Limites_Sectores_N(j+1)))) / Muestras_Sectores_N(j);
end

% k_j

k_j = zeros(1, N);
for j = 1:N
    
    k_j(j) = (23.29041409 - 16.8617370 * sqrt(s_j(j)^2 + c_j(j)^2) - 17.4749884 * exp(-(s_j(j)^2 + c_j(j)^2)))^(-1);
end

% mu_j 
mu_j = zeros(1, N);
for j = 1:N
    
    s_j_val = s_j(j);
    c_j_val = c_j(j);
    if  s_j_val>=0 && c_j_val > 0
        mu_j(j) = atan(s_j_val / c_j_val);

    elseif s_j_val > 0 && c_j_val == 0 
        mu_j(j) = pi / 2;

    elseif  c_j_val <= 0 
        mu_j(j) = pi +atan(s_j_val / c_j_val);

    elseif s_j_val == 0 && c_j_val == -1 
        mu_j(j) = pi;

    elseif s_j_val < 0 && c_j_val > 0 
        mu_j(j) = 2*pi +atan(s_j_val / c_j_val);

    elseif s_j_val < 0 && c_j_val == 0 
       mu_j(j) = 3*pi/2;       

    end
end


x0 = [zeros(1, N), k_j,mu_j];
lb = [zeros(1, N), zeros(1, N), zeros(1, N)];
ub = [ones(1, N), Inf * ones(1, N),2*pi*ones(1, N)];


Aeq = [ones(1, N),zeros(1, 2*N)];
beq = 1;

% Opciones para lsqnonlin
options = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt','Display', 'iter');

x = lsqnonlin(@(x)(S(P,T,N,Limites_Sectores_T,x)), x0, lb, ub,[],[],Aeq, beq, options);


function y = S(P, T, N, Limites_Sectores_T, x)
    y = 0;

    for i = 1:T-1
        Z = 0;

        for j = 1:N
            Z = Z + x(j) * integral(@(theta) exp((x(j+N) * cos(theta - x(j+2*N))) / (2*pi*besseli(0, x(j+N)))), 0, Limites_Sectores_T(i));
        end

        y = y + (P(i) - Z)^2;
    end
end

This is the error message I get:

Error using Distribucion_WD>@(x)(S(P,T,N,Limites_Sectores_T,x))
Too many input arguments.

Error in lsqnonlin (line 218)
            initVals.F = feval(funfcn{3},xCurrent,varargin{:});

Error in Distribucion_WD (line 110)
x = lsqnonlin(@(x)(S(P,T,N,Limites_Sectores_T,x)), x0, lb, ub,[],[],Aeq, beq, options);

Caused by:
    Failure in initial objective function evaluation. LSQNONLIN cannot continue.

What am I doing wrong? It looks like I'm passing too many arguments, but I'm passing the same arguments that I defined in the function.

I have tried the following

x = lsqnonlin(@(x)(S(P, N, Limites_Sectores_T, x)), x0, lb, ub, Aeq, beq,optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt','Display', 'iter'));

In this case, i have the following error.

     Error using lsqnonlin
Invalid datatype. Options argument must be created with OPTIMOPTIONS.

Error in Distribucion_WD (line 112)
 x = lsqnonlin(@(x)(S(P, N, Limites_Sectores_T, x)), x0, lb, ub, Aeq, beq,optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt','Display', 'iter'));

In short, how should I define the objective function to avoid all these problems? Why do I get these errors? Thank you very much for your time