Matlab Genetic Algorithm for a non-genetic case

Question

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!

It sounds to me like you're trying to find the least squares solution for a linear system of equations. If so, you don't need to explicitly optimize, because there are analytic solutions for such problems. Basically, that's what the field of linear algebra is about. – Whether my interpretation is correct or not, the details of your question are pretty unclear. Can your provide some code *that is self-contained so we can run it by ourselves*, and then explain your problem with respect to that code? — A. Donda, Apr 06 '15 at 14:22
@A. Donda - As you asked I tried to explain better my question, thank you very much for your answer — PingP, Apr 06 '15 at 15:37
Thanks, but I ran your code several times, and the RMSE is always 0 – it cannot get any better than this. Doesn't this mean you have achieved what you wanted? — A. Donda, Apr 06 '15 at 15:42
@A. Donda - I had to update the code, I forgot to put a value in a field of the 'randi' — PingP, Apr 06 '15 at 15:44
I'm still getting 0 all the time. Also, why do you have so many `randi`s in your code? Are you simulating data, or is this part of your attempt to optimize? — A. Donda, Apr 06 '15 at 15:47
OK, so you have `Gparque` that depends on `pares`, which are 27 x 2 variables from the range [0, 90], and you want to find the value of `pares` such that `Gparque` matches veic. Correct? Can the values of `pares` be only integers? — A. Donda, Apr 06 '15 at 15:53
Yes, I was just about to say that in my code the RMSE is 0,43; the randomness is just to simulate data — PingP, Apr 06 '15 at 15:54
If the `pares` can only be integers, that's a problem, because most optimization algorithms are designed for continuous spaces. — A. Donda, Apr 06 '15 at 16:16
@A. Donda - They can be integers, but it is not mandatory that they are — PingP, Apr 06 '15 at 16:19
Ok, and what about the constraint to [0, 90]. Can that be lifted, too? — A. Donda, Apr 06 '15 at 16:34
Is it ok to have values in pares smaller than 0 or larger than 90? — A. Donda, Apr 06 '15 at 16:48
It is ok but it isn't ideal, but do you know anyway to resolve this doubt? Even if it means change the range of the values of pares? — PingP, Apr 06 '15 at 16:53

A. Donda · Accepted Answer · 2015-04-06T17:46:02.087

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.

I've tried and I think it is working. But still some problems: 1. I wanted more iterations and I was trying: options = optimset('MaxFunEvals',10000); pares_opt = fminsearch(@FUN, randi([0,90],27,2),options), but it isn't working, do you know another way? 2. Sometimes the evaluations are not enough, so when I dont need more iterations it doesn't return the optimum RMSE. How can I see these value, the RMSE that corresponds to the values of pares that are returned? 3. And just one more question, pares_opt changes values everytime I run, any suggestion? Thank you so much again! — PingP, Apr 06 '15 at 22:09
@PingP, 1. Try setting both `MaxFunEvals` and `MaxIter`. 2. Just compute `RMSE = fun(pares_opt)`. 3. Yes, I mentioned that in the answer. I think it is because for those random simulated data, there is no clearly optimal solution, but lots of local minima. Try running the code on real data. — A. Donda, Apr 08 '15 at 02:01
I have that problem with the real data, that's my main issue and I have already set MaxFunEvals and MaxIter but it still sometimes says that it is not enough and sometimes it is.. The random matrix is the problem I think, I need to modify that part, that's why I think that genetic algorithm maybe the best way to approach this — PingP, Apr 08 '15 at 08:57
@PingP I see. Something else you can try is to run the algorithm many times with different random initialization, save the solution and associated RMSE, and in the end pick the version with the smallest RMSE. Not perfect, but you'll get a better chance at finding the best local minimum. Beyond that I can't really help you any more because I do not understand the underlying problem. I guess it has something to do with astronomy? — A. Donda, Apr 08 '15 at 14:03

Matlab Genetic Algorithm for a non-genetic case

1 Answers1