# How to optimise an objective function with a summation of integrals

조회 수: 4 (최근 30일)
Samuel M 2024년 3월 15일
편집: Torsten 2024년 3월 19일
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_clean = WD(~isnan(WD));
T= 360;
P = cumsum(Muestras_Sectores_T) / total_samples;
N = 4;
% s_j y c_j
s_j = zeros(1, N);
c_j = zeros(1, N);
for j = 1:N
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 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 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 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이전 댓글 -1개 표시이전 댓글 -1개 숨기기
Samuel M 2024년 3월 18일
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 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.
##### 댓글 수: 2없음 표시없음 숨기기
Samuel M 2024년 3월 19일
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 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 CenterFile Exchange에서 Numerical Integration and Differential Equations에 대해 자세히 알아보기

R2022a

### Community Treasure Hunt

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

Start Hunting!

Translated by