How to optimise an objective function with a summation of integrals

2 次查看（过去 30 天）
Samuel M
Samuel M 2024-3-15
编辑： Torsten 2024-3-19
Ran in:
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
  2 个评论
Torsten
Torsten 2024-3-15
编辑：Torsten 2024-3-15
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
Samuel M 2024-3-18
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

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Walter Roberson
Walter Roberson 2024-3-16
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
  1 个评论
Samuel M
Samuel M 2024-3-18
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.

请先登录，再进行评论。

Torsten
Torsten 2024-3-18
移动：Walter Roberson 2024-3-18
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.
  2 个评论
Samuel M
Samuel M 2024-3-19
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
Torsten 2024-3-19
编辑：Torsten 2024-3-19
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Get Started with Optimization Toolbox 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by