How to optimise an objective function with a summation of integrals
조회 수: 3 (최근 30일)
이전 댓글 표시
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);
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.
답변 (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
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.
참고 항목
카테고리
Help Center 및 File Exchange에서 Numerical Integration and Differential Equations에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!