scilab/modules/genetic_algorithms/macros/optim_nsga.sci

   1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   2 // Copyright (C) 2008 - Yann COLLETTE <yann.collette@renault.com>
   3 //
   4 // This file must be used under the terms of the CeCILL.
   5 // This source file is licensed as described in the file COPYING, which
   6 // you should have received as part of this distribution.  The terms
   7 // are also available at
   8 // http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
   9
  10 function [pop_opt, fobj_pop_opt, pop_init, fobj_pop_init] = optim_nsga(ga_f, pop_size, nb_generation, p_mut, p_cross, Log, param, sigma, pow)
  11
  12     [nargout, nargin] = argn();
  13
  14     if ~isdef("param","local") then
  15         param = [];
  16     end
  17
  18     [codage_func,err]    = get_param(param,"codage_func",coding_ga_identity);
  19     [init_func,err]      = get_param(param,"init_func",init_ga_default);
  20     [crossover_func,err] = get_param(param,"crossover_func",crossover_ga_default);
  21     [mutation_func,err]  = get_param(param,"mutation_func",mutation_ga_default);
  22     [selection_func,err] = get_param(param,"selection_func",selection_ga_elitist);
  23     [nb_couples,err]     = get_param(param,"nb_couples",100);
  24     [pressure,err]       = get_param(param,"pressure",0.05);
  25
  26     if ~isdef("ga_f","local") then
  27         error(gettext("optim_moga: ga_f is mandatory"));
  28     else
  29         if typeof(ga_f)=="list" then
  30             deff("y=_ga_f(x)","y=ga_f(1)(x, ga_f(2:$))");
  31         else
  32             deff("y=_ga_f(x)","y=ga_f(x)");
  33         end
  34     end
  35
  36     if ~isdef("pop_size","local") then
  37         pop_size = 100;
  38     end
  39     if ~isdef("nb_generation","local") then
  40         nb_generation = 10;
  41     end
  42     if ~isdef("p_mut","local") then
  43         p_mut = 0.1;
  44     end
  45     if ~isdef("p_cross","local") then
  46         p_cross = 0.1;
  47     end
  48     if ~isdef("Log","local") then
  49         Log = %F;
  50     end
  51     if ~isdef("sigma","local") then
  52         sigma = 0.01;
  53     end
  54     if ~isdef("pow","local") then
  55         pow = 2;
  56     end
  57
  58     // Initialization of the population
  59     if (Log) then
  60         printf(gettext("%s: Initialization of the population\n"),"optim_nsga");
  61     end
  62
  63     Pop = list();
  64     Pop = init_func(pop_size,param);
  65
  66
  67     if (nargout>=3) then
  68         pop_init = Pop;
  69     end
  70
  71     // Code the individuals
  72     Pop = codage_func(Pop,"code",param);
  73
  74     for i=1:length(Pop)
  75         MO_FObj_Pop(i,:) = _ga_f(Pop(i));
  76     end
  77
  78     // Compute the domination rank
  79     for i=1:size(MO_FObj_Pop,1)
  80         Index = 0;
  81         for j=1:size(MO_FObj_Pop,1)
  82             Index = Index + double(and(MO_FObj_Pop(i,:)<=MO_FObj_Pop(j,:)) & or(MO_FObj_Pop(i,:)<MO_FObj_Pop(j,:)));
  83         end
  84         FObj_Pop(i) = - (Index + 1);
  85     end
  86
  87     if (nargout==4) then
  88         fobj_pop_init = MO_FObj_Pop;
  89     end
  90
  91     // The genetic algorithm
  92     for i=1:nb_generation
  93         if (Log) then
  94             printf(gettext("%s: iteration %d / %d"), "optim_nsga", i, nb_generation);
  95         end
  96
  97         // Computation of the niching penality
  98         for j=1:size(MO_FObj_Pop,1)
  99             Niching(j) = 0;
 100             for k=1:size(MO_FObj_Pop,1)
 101                 Distance = sqrt(sum((MO_FObj_Pop(j,:) - MO_FObj_Pop(k,:)).^2));
 102                 if Distance < sigma then
 103                     Niching(j) = Niching(j) + (1 - Distance / sigma)^pow;
 104                 end
 105             end
 106         end
 107
 108         // Apply niching penality to fobj
 109         FObj_Pop = FObj_Pop ./ Niching;
 110
 111         FObj_Pop_Max = max(FObj_Pop);
 112         FObj_Pop_Min = min(FObj_Pop);
 113
 114         // Normalization of the efficiency
 115         Efficiency = (1 - pressure) * (FObj_Pop_Max - FObj_Pop) / max([FObj_Pop_Max - FObj_Pop_Min %eps]) + pressure;
 116
 117         //
 118         // Selection
 119         //
 120         Indiv1 = list();
 121         Indiv2 = list();
 122         Wheel = cumsum(Efficiency);
 123
 124         for j=1:nb_couples
 125             // Selection of the first individual in the couple
 126             Shoot = grand(1,1,"def")*Wheel($);
 127             Index = find(Wheel<=Shoot);
 128             if length(Index)>1 then Index = Index($); end;
 129             if isempty(Index)  then Index = 1; end;
 130             Indiv1(j)           = Pop(Index);
 131             MO_FObj_Indiv1(j,:) = MO_FObj_Pop(Index,:);
 132             // Selection of the second individual in the couple
 133             Shoot = grand(1,1,"def")*Wheel($);
 134             Index = find(Wheel<=Shoot);
 135             if length(Index)>1 then Index = Index($); end;
 136             if isempty(Index)  then Index = 1; end;
 137             Indiv2(j)           = Pop(Index);
 138             MO_FObj_Indiv2(j,:) = MO_FObj_Pop(Index,:);
 139         end
 140         //
 141         // Crossover
 142         //
 143         for j=1:nb_couples
 144             if (p_cross>grand(1,1,"def")) then
 145                 [x1, x2]  = crossover_func(Indiv1(j), Indiv2(j),param);
 146                 Indiv1(j) = x1;
 147                 Indiv2(j) = x2;
 148                 ToCompute_I1(j) = %T;
 149                 ToCompute_I2(j) = %T;
 150             else
 151                 ToCompute_I1(j) = %F;
 152                 ToCompute_I2(j) = %F;
 153             end
 154         end
 155         //
 156         // Mutation
 157         //
 158         for j=1:nb_couples
 159             if (p_mut>grand(1,1,"def")) then
 160                 x1 = mutation_func(Indiv1(j),param);
 161                 Indiv1(j) = x1;
 162                 ToCompute_I1(j) = %T;
 163             end
 164             if (p_mut>grand(1,1,"def")) then
 165                 x2 = mutation_func(Indiv2(j),param);
 166                 Indiv2(j) = x2;
 167                 ToCompute_I2(j) = %T;
 168             end
 169         end
 170         //
 171         // Computation of the objective functions
 172         //
 173         for j=1:length(Indiv1)
 174             if ToCompute_I1(j) then MO_FObj_Indiv1(j,:) = _ga_f(Indiv1(j)); end
 175             if ToCompute_I2(j) then MO_FObj_Indiv2(j,:) = _ga_f(Indiv2(j)); end
 176         end
 177
 178         // Reinit ToCompute lists
 179         ToCompute_I1 = ToCompute_I1 & %F;
 180         ToCompute_I2 = ToCompute_I2 & %F;
 181
 182         // Compute the domination rank
 183         for j=1:size(MO_FObj_Indiv1,1)
 184             // We compute the rank for Indiv1
 185             Index1 = 0; Index2 = 0; Index3 = 0;
 186             for k=1:size(MO_FObj_Indiv1,1)
 187                 Index1 = Index1 + double(and(MO_FObj_Indiv1(j,:)<=MO_FObj_Indiv1(k,:)) & or(MO_FObj_Indiv1(j,:)<MO_FObj_Indiv1(k,:)));
 188                 Index2 = Index2 + double(and(MO_FObj_Indiv1(j,:)<=MO_FObj_Indiv2(k,:)) & or(MO_FObj_Indiv1(j,:)<MO_FObj_Indiv2(k,:)));
 189             end
 190             for k=1:size(MO_FObj_Pop,1)
 191                 Index3 = Index3 + and(MO_FObj_Indiv1(j,:)<=MO_FObj_Pop(k,:)) & or(MO_FObj_Indiv1(j,:)<MO_FObj_Pop(k,:));
 192             end
 193             FObj_Indiv1(j) = - (Index1 + Index2 + Index3 + 1);
 194
 195             // We compute the rank for Indiv2
 196             Index1 = 0; Index2 = 0; Index3 = 0;
 197             for k=1:size(MO_FObj_Indiv1,1)
 198                 Index1 = Index1 + double(and(MO_FObj_Indiv2(j,:)<=MO_FObj_Indiv1(k,:)) & or(MO_FObj_Indiv2(j,:)<MO_FObj_Indiv1(k,:)));
 199                 Index2 = Index2 + double(and(MO_FObj_Indiv2(j,:)<=MO_FObj_Indiv2(k,:)) & or(MO_FObj_Indiv2(j,:)<MO_FObj_Indiv2(k,:)));
 200             end
 201             for k=1:size(MO_FObj_Pop,1)
 202                 Index3 = Index3 + double(and(MO_FObj_Indiv2(j,:)<=MO_FObj_Pop(k,:)) & or(MO_FObj_Indiv2(j,:)<MO_FObj_Pop(k,:)));
 203             end
 204             FObj_Indiv2(j) = - (Index1 + Index2 + Index3 + 1);
 205         end
 206
 207         // We compute the rank for Pop
 208         for j=1:size(MO_FObj_Pop,1)
 209             Index1 = 0; Index2 = 0; Index3 = 0;
 210             for k=1:size(MO_FObj_Indiv1,1)
 211                 Index1 = Index1 + double(and(MO_FObj_Pop(j,:)<=MO_FObj_Indiv1(k,:)) & or(MO_FObj_Pop(j,:)<MO_FObj_Indiv1(k,:)));
 212                 Index2 = Index2 + double(and(MO_FObj_Pop(j,:)<=MO_FObj_Indiv2(k,:)) & or(MO_FObj_Pop(j,:)<MO_FObj_Indiv2(k,:)));
 213             end
 214             for k=1:size(FObj_Pop,1)
 215                 Index3 = Index3 + double(and(MO_FObj_Pop(j,:)<=MO_FObj_Pop(k,:)) & or(MO_FObj_Pop(j,:)<MO_FObj_Pop(k,:)));
 216             end
 217             FObj_Pop(j) = - (Index1 + Index2 + Index3 + 1);
 218         end
 219
 220         //
 221         // Recombination
 222         //
 223
 224         [Pop,FObj_Pop,Efficiency,MO_FObj_Pop] = selection_func(Pop,Indiv1,Indiv2,FObj_Pop,FObj_Indiv1,FObj_Indiv2, ...
 225         MO_FObj_Pop,MO_FObj_Indiv1,MO_FObj_Indiv2,param);
 226         if (Log) then
 227             printf(gettext(" - min / max value found = %f / %f\n"), min(FObj_Pop), max(FObj_Pop));
 228         end
 229     end
 230
 231     pop_opt      = codage_func(Pop,"decode",param);
 232     fobj_pop_opt = MO_FObj_Pop;
 233 endfunction