How to optimise an objective function with a summation of integrals

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:

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

filenames = unzip('dates.zip');
WD = ncread(filenames{1,1}, 'WD');


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

T= 360;

Limites_Sectores_T = deg2rad(0:360/T:360); 


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 = 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 = [0.5*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);

Warning: Constraints not supported by selected algorithm. Switching algorithm to interior-point.

Initial point X0 is not between bounds LB and UB; 
FMINCON shifted X0 to strictly satisfy the bounds.

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0      13    3.470750e+06    8.000e-01    7.621e+04
    1      26    8.943390e+05    0.000e+00    1.149e+06    9.791e+01
    2      39    8.940652e+05    0.000e+00    1.031e+06    2.128e+00
    3      53    8.743310e+05    0.000e+00    9.369e+05    2.722e+00
    4      68    8.728424e+05    0.000e+00    9.341e+05    1.505e+00
    5      81    8.492274e+05    0.000e+00    4.951e+05    1.672e+00
    6      96    8.456461e+05    0.000e+00    4.875e+05    8.171e-01
    7     110    8.426631e+05    0.000e+00    4.808e+05    1.163e+00
    8     123    8.363319e+05    0.000e+00    4.223e+05    2.790e+00
    9     136    8.297205e+05    0.000e+00    3.939e+05    7.468e-01
   10     149    8.288938e+05    0.000e+00    3.640e+05    1.109e-01
   11     162    8.259818e+05    0.000e+00    2.620e+04    5.843e-01
   12     175    8.258132e+05    0.000e+00    6.283e+03    3.625e-02
   13     188    8.257974e+05    0.000e+00    1.372e+03    2.546e-02
   14     201    8.257929e+05    1.110e-16    1.294e+01    4.631e-03
   15     214    8.257929e+05    0.000e+00    4.672e-01    6.963e-05

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

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);
    end
end

However in Matlab version 2022a this is the error I get

Error using Distribucion_WD>@(x)(S(P,T,N,Limites_Sectores_T,x))

Too many input arguments.

What am I doing wrong?

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

[Torsten]

This does not explain the error, but if you use "lsqnonlin", you have to return the T terms

P_k - sum_j(...) (1 <= k <= T)

separately, not the sum of squares

sum_k (P_k - sum_j(...))^2

in one single value y.

You should provide executable code that reproduces the error message you get. Above, at least the .nc file is missing.

[Samuel M]

Hello Torsten.

Thank you for your help. I have already changed the function definition as you indicated.

I have added the data file for the check as you suggested.

If I run the code here on the forum with version 2023. The message I get is the one you see now. However, in my version of Matlab 2022, I still get the error Too many arguments

Answer 1

[Walter Roberson]

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

There are a few different valid calling forms for lsqnonlin(), but if you go as far as Aeq beq then the next parameter must be nonlcon before options.

You can only short-circuit options if you place it directly after lb, ub

[Samuel M]

Thank you for your reply Walter.

I have modified what you said.

But I still seem to have the same problem.

I have uploaded the data file to the discussion and the code appearance is now as shown.

With Matlab 2023b, I don't seem to have the problem I describe, but with Matlab 2022a, I still get the Too many arguments error.

Answer 2

[Torsten]

The support of Aeq and beq inputs was introduced in release R2023a.

You should always consult the documentation relevant for your release, not the most recent one.

I guess you will have to switch to "fmincon" if you want to keep the constraint.

R2023a: Linear and Nonlinear Constraint Support

lsqnonlin gains support for both linear and nonlinear constraints. To enable constraint satisfaction, the solver uses the "interior-point" algorithm from fmincon.

• If you specify constraints but do not specify an algorithm, the solver automatically switches to the "interior-point" algorithm.
• If you specify constraints and an algorithm, you must specify the "interior-point" algorithm.

For algorithm details, see Modified fmincon Algorithm for Constrained Least Squares. For an example, see Compare lsqnonlin and fmincon for Constrained Nonlinear Least Squares.

[Samuel M]

Thank you for your response.

I will try to solve it with fmincon. I was only using lsqnonlin because it has the option to use the Levenberg-Marquardt algorithm, which is the one used in the paper.

[Torsten]

Note that in "fmincon", you have to use your old formulation of the objective function, thus

y = sum_k (P_k - sum_j(...))^2

in one single value y.