I've never used optimization tools, but I think I have to use now, so I'm a bit lost. After using the answer given by @A. Donda, I have noticed that maybe that is not the best solution because every time I run the function it gives a different matrix 'pares' and in the majority of times it says that I need more evaluations. So I was thinking that maybe the answer to my problem are Genetic Algorithms optimization, but once again I do not know how to work with them.
My first problem is described below and the answer by @A. Donda is in the only post of a answer. I really need this optimization done and I don't know how to proceed for this case with GA tools.
Thank you so much in advance again, and thank you @A. Donda for your answer.
As asked, I tried to put here the code that I was trying to explain, I hope it will result:
function opt_pares
clear all; clc; close all;
h = randi(24,8760,1);
nd = randi(365,8760,1);
veic = randi(333,8760,1);
max_veic = max(veic);
veicN = veic./max_veic;
Gh = randi(500,8760,1);
Dh = randi(500,8760,1);
Ih = Gh-Dh;
A = randi([300 800], 27,1);
max_Gh = max(Gh);
max_Dh = max(Dh);
max_Ih = max(Ih);
lat = 70;
HRA =15.*(h-12);
decl = 23.27*sind(360*(284+nd)/365);
Ii = zeros(8760,27);
Di = zeros(8760,27);
Gi = zeros(8760,27);
pares = randi([0,90],27,2);
inclin = pares(:,1);
azim = pares(:,2);
% for MATRIZC
for n=1:27
Ii(:,n) = Ih.*(sind(decl).*sind(lat).*cosd(inclin(n))-sind(decl).*cosd(lat).*sind(inclin(n)).*cosd(azim(n))+cosd(decl).*cosd(lat).*cosd(inclin(n)).*cosd(HRA)+cosd(decl).*sind(lat).*sind(inclin(n)).*cosd(azim(n)).*cosd(HRA)+cosd(decl).*sind(inclin(n)).*sind(azim(n)).*sind(HRA));
Di(:,n) = 0.5*Dh.*(1+cosd(inclin(n)));
Gi(:,n) = (Ii(:,n)+Di(:,n))*A(n,1);
end
Gparque = sum(Gi,2);
max_Gparque = max(Gparque);
GparqueN = Gparque./max_Gparque;
RMSE = sqrt(mean((GparqueN-veicN).^2));
% end
end
I don't know if it is possible, maybe this time I can be more assertive.
My main goal is to achieve the best 'RMSE' possible, to do so I have to create a matrix ('pares') where each line contains a pair of values (one value from each column).
These values have to be within a certain range(0-90). With each of this 27 pairs I have to calculate 'Ii'/'Gi'/'Di', giving me a matrix with a size like 8760*27.
Then I make a sum of 'Gi' to have 'Gparque'(vector 8760*1) and finally I I calculate 'RMSE'. When I have RMSE calculated, I have to modify the matrix 'pares' to other values that can result in a better RMSE. Once there are many combinations of 27 values that can be within the 0-90 range, I have to get a solution that can optimize this search for the minimum RMSE.
The parts that are in comments in the code (a for loop with 'pares') is the thing that I have no idea how to do, because I have to change the values of 'pares' but with some optimization criteria that can approximate the minimum of RMSE.
I hope this time I have explain this doubt better.
Thank you very much!
OK, so here is an attempt at a question. I'm not sure how useful the results will be in the end, because I don't understand the underlying problem and I don't have real data to test it with.
You were right that you need an optimization algorithm, your problem appears to be more complex than simple linear algebra. For optimization I use the function fminsearch from the Optmization Toolbox.
First the function whose value is to be optimized (the objective function) needs to be defined. Based on your code, this is
function RMSE = fun(pares)
inclin = pares(:,1);
azim = pares(:,2);
Ii = zeros(8760,27);
Di = zeros(8760,27);
Gi = zeros(8760,27);
for n=1:27
Ii(:,n) = Ih.*(sind(decl).*sind(lat).*cosd(inclin(n))-sind(decl).*cosd(lat).*sind(inclin(n)).*cosd(azim(n))+cosd(decl).*cosd(lat).*cosd(inclin(n)).*cosd(HRA)+cosd(decl).*sind(lat).*sind(inclin(n)).*cosd(azim(n)).*cosd(HRA)+cosd(decl).*sind(inclin(n)).*sind(azim(n)).*sind(HRA));
Di(:,n) = 0.5*Dh.*(1+cosd(inclin(n)));
Gi(:,n) = (Ii(:,n)+Di(:,n))*A(n,1);
end
Gparque = sum(Gi,2);
max_Gparque = max(Gparque);
GparqueN = Gparque./max_Gparque;
RMSE = sqrt(mean((GparqueN-veicN).^2));
end
Now we can call
pares_opt = fminsearch(@fun, randi([0,90],27,2))
using random initialization. The optimization takes quite a while because the objective function is not very efficiently implemented. Here's a vectorized version that does the same:
% precompute
cHRA = cosd(HRA);
sHRA = sind(HRA);
sdecl = sind(decl);
cdecl = cosd(decl);
slat = sind(lat);
clat = cosd(lat);
function RMSE = fun(pares)
% precompute
cinclin = cosd(pares(:,1))';
sinclin = sind(pares(:,1))';
cazim = cosd(pares(:,2))';
sazim = sind(pares(:,2))';
Ii = bsxfun(@times, Ih, ...
sdecl * (slat * cinclin - clat * sinclin .* cazim) ...
+ (cdecl .* cHRA) * (clat * cinclin + slat * sinclin .* cazim) ...
+ (cdecl .* sHRA) * (sinclin .* sazim));
Di = 0.5 * Dh * (1 + cinclin);
Gi = (Ii + Di) * diag(A);
Gparque = sum(Gi,2);
max_Gparque = max(Gparque);
GparqueN = Gparque./max_Gparque;
RMSE = sqrt(mean((GparqueN-veicN).^2));
end
We have not yet implemented the constraint for pares to lie within [0, 90]. A crude way to do this is to insert these lines:
if any(pares(:) < 0) || any(pares(:) > 90)
RMSE = inf;
return
end
at the beginning of the objective function.
Putting it all together:
function Raquel
h = randi(24,8760,1);
nd = randi(365,8760,1);
veic = randi(333,8760,1);
max_veic = max(veic);
veicN = veic./max_veic;
Gh = randi(500,8760,1);
Dh = randi(500,8760,1);
Ih = Gh-Dh;
A = randi([300 800], 27,1);
lat = 70;
HRA =15.*(h-12);
decl = 23.27*sind(360*(284+nd)/365);
% precompute
cHRA = cosd(HRA);
sHRA = sind(HRA);
sdecl = sind(decl);
cdecl = cosd(decl);
slat = sind(lat);
clat = cosd(lat);
pares_opt = fminsearch(@fun, randi([0,90],27,2))
function RMSE = fun(pares)
if any(pares(:) < 0) || any(pares(:) > 90)
RMSE = inf;
return
end
% precompute
cinclin = cosd(pares(:,1))';
sinclin = sind(pares(:,1))';
cazim = cosd(pares(:,2))';
sazim = sind(pares(:,2))';
Ii = bsxfun(@times, Ih, ...
sdecl * (slat * cinclin - clat * sinclin .* cazim) ...
+ (cdecl .* cHRA) * (clat * cinclin + slat * sinclin .* cazim) ...
+ (cdecl .* sHRA) * (sinclin .* sazim));
Di = 0.5 * Dh * (1 + cinclin);
Gi = (Ii + Di) * diag(A);
Gparque = sum(Gi,2);
max_Gparque = max(Gparque);
GparqueN = Gparque./max_Gparque;
RMSE = sqrt(mean((GparqueN-veicN).^2));
end
end
With simulated data, if I run the optimization twice on the same randomized data but different initial values I get different solutions. This is an indication that there is more than one local minimum of the objective function. Hopefully, this will not be the case with real data.