1  using System;


2  using System.Collections.Generic;


3  using System.Diagnostics;


4  using System.Linq;


5  using System.Threading;


6  using HEAL.Attic;


7  using HeuristicLab.Common;


8  using HeuristicLab.Core;


9  using HeuristicLab.Data;


10  using HeuristicLab.Encodings.RealVectorEncoding;


11  using HeuristicLab.Optimization;


12  using HeuristicLab.Parameters;


13  using HeuristicLab.Random;


14 


15  namespace HeuristicLab.Algorithms.NSGA3


16  {


17  /// <summary>


18  /// The Reference Point Based Nondominated Sorting Genetic Algorithm III was introduced in Deb


19  /// et al. 2013. An Evolutionary ManyObjective Optimization Algorithm Using Reference Point


20  /// Based Nondominated Sorting Approach. IEEE Transactions on Evolutionary Computation, 18(4),


21  /// pp. 577601.


22  /// </summary>


23  [Item("NSGAIII", "The Reference Point Based Nondominated Sorting Genetic Algorithm III was introduced in Deb et al. 2013. An Evolutionary ManyObjective Optimization Algorithm Using Reference Point Based Nondominated Sorting Approach. IEEE Transactions on Evolutionary Computation, 18(4), pp. 577601.")]


24  [Creatable(Category = CreatableAttribute.Categories.PopulationBasedAlgorithms, Priority = 136)]


25  [StorableType("07C745F7A8A34F998B2CF97E639F9AC3")]


26  public class NSGA3 : BasicAlgorithm


27  {


28  private const double EPSILON = 10e6; // a tiny number that is greater than 0


29 


30  public override bool SupportsPause => false;


31 


32  #region ProblemProperties


33 


34  public override Type ProblemType


35  {


36  get { return typeof(MultiObjectiveBasicProblem<RealVectorEncoding>); }


37  }


38 


39  public new MultiObjectiveBasicProblem<RealVectorEncoding> Problem


40  {


41  get { return (MultiObjectiveBasicProblem<RealVectorEncoding>)base.Problem; }


42  set { base.Problem = value; }


43  }


44 


45  public int NumberOfObjectives => Problem.Maximization.Length;


46 


47  #endregion ProblemProperties


48 


49  #region Storable fields


50 


51  [Storable]


52  private IRandom random;


53 


54  [Storable]


55  private List<Solution> solutions;


56 


57  #endregion Storable fields


58 


59  #region ParameterAndResultsNames


60 


61  // Parameter Names


62 


63  private const string SeedName = "Seed";


64  private const string SetSeedRandomlyName = "SetSeedRandomly";


65  private const string PopulationSizeName = "PopulationSize";


66  private const string CrossoverProbabilityName = "CrossoverProbability";


67  private const string CrossoverContiguityName = "CrossoverContiguity";


68  private const string MutationProbabilityName = "MutationProbability";


69  private const string MaximumGenerationsName = "MaximumGenerations";


70  private const string DominateOnEqualQualitiesName = "DominateOnEqualQualities";


71 


72  // Results Names


73 


74  private const string GeneratedReferencePointsResultName = "Generated Reference Points";


75  private const string CurrentGenerationResultName = "Generations";


76  private const string CurrentFrontResultName = "Pareto Front"; // Do not touch this


77 


78  #endregion ParameterAndResultsNames


79 


80  #region ParameterProperties


81 


82  private IFixedValueParameter<IntValue> SeedParameter


83  {


84  get { return (IFixedValueParameter<IntValue>)Parameters[SeedName]; }


85  }


86 


87  private IFixedValueParameter<BoolValue> SetSeedRandomlyParameter


88  {


89  get { return (IFixedValueParameter<BoolValue>)Parameters[SetSeedRandomlyName]; }


90  }


91 


92  private IFixedValueParameter<IntValue> PopulationSizeParameter


93  {


94  get { return (IFixedValueParameter<IntValue>)Parameters[PopulationSizeName]; }


95  }


96 


97  private IFixedValueParameter<PercentValue> CrossoverProbabilityParameter


98  {


99  get { return (IFixedValueParameter<PercentValue>)Parameters[CrossoverProbabilityName]; }


100  }


101 


102  private IFixedValueParameter<DoubleValue> CrossoverContiguityParameter


103  {


104  get { return (IFixedValueParameter<DoubleValue>)Parameters[CrossoverContiguityName]; }


105  }


106 


107  private IFixedValueParameter<PercentValue> MutationProbabilityParameter


108  {


109  get { return (IFixedValueParameter<PercentValue>)Parameters[MutationProbabilityName]; }


110  }


111 


112  private IFixedValueParameter<IntValue> MaximumGenerationsParameter


113  {


114  get { return (IFixedValueParameter<IntValue>)Parameters[MaximumGenerationsName]; }


115  }


116 


117  private IFixedValueParameter<BoolValue> DominateOnEqualQualitiesParameter


118  {


119  get { return (IFixedValueParameter<BoolValue>)Parameters[DominateOnEqualQualitiesName]; }


120  }


121 


122  #endregion ParameterProperties


123 


124  #region Properties


125 


126  public IntValue Seed => SeedParameter.Value;


127 


128  public BoolValue SetSeedRandomly => SetSeedRandomlyParameter.Value;


129 


130  public IntValue PopulationSize => PopulationSizeParameter.Value;


131 


132  public PercentValue CrossoverProbability => CrossoverProbabilityParameter.Value;


133 


134  public DoubleValue CrossoverContiguity => CrossoverContiguityParameter.Value;


135 


136  public PercentValue MutationProbability => MutationProbabilityParameter.Value;


137 


138  public IntValue MaximumGenerations => MaximumGenerationsParameter.Value;


139 


140  public BoolValue DominateOnEqualQualities => DominateOnEqualQualitiesParameter.Value;


141 


142  public List<List<Solution>> Fronts { get; private set; }


143 


144  public List<ReferencePoint> ReferencePoints { get; private set; }


145 


146  // todo: create one property for the Generated Reference Points and one for the current


147  // generations reference points


148 


149  #endregion Properties


150 


151  #region ResultsProperties


152 


153  public DoubleMatrix ResultsGeneratedReferencePoints


154  {


155  get { return (DoubleMatrix)Results[GeneratedReferencePointsResultName].Value; }


156  set { Results[GeneratedReferencePointsResultName].Value = value; }


157  }


158 


159  public DoubleMatrix ResultsSolutions


160  {


161  get { return (DoubleMatrix)Results[CurrentFrontResultName].Value; }


162  set { Results[CurrentFrontResultName].Value = value; }


163  }


164 


165  public IntValue ResultsCurrentGeneration


166  {


167  get { return (IntValue)Results[CurrentGenerationResultName].Value; }


168  set { Results[CurrentGenerationResultName].Value = value; }


169  }


170 


171  #endregion ResultsProperties


172 


173  #region Constructors


174 


175  public NSGA3() : base()


176  {


177  Parameters.Add(new FixedValueParameter<IntValue>(SeedName, "The random seed used to initialize the new pseudo random number generator.", new IntValue(0)));


178  Parameters.Add(new FixedValueParameter<BoolValue>(SetSeedRandomlyName, "True if the random seed should be set to a random value, otherwise false.", new BoolValue(true)));


179  Parameters.Add(new FixedValueParameter<IntValue>(PopulationSizeName, "The size of the population of Individuals.", new IntValue(200)));


180  Parameters.Add(new FixedValueParameter<PercentValue>(CrossoverProbabilityName, "The probability that the crossover operator is applied on two parents.", new PercentValue(0.9)));


181  Parameters.Add(new FixedValueParameter<DoubleValue>(CrossoverContiguityName, "The contiguity value for the Simulated Binary Crossover that specifies how close a child should be to its parents (larger value means closer). The value must be greater than or equal than 0. Typical values are in the range [2;5]."));


182  Parameters.Add(new FixedValueParameter<PercentValue>(MutationProbabilityName, "The probability that the mutation operator is applied on a Individual.", new PercentValue(0.05)));


183  Parameters.Add(new FixedValueParameter<IntValue>(MaximumGenerationsName, "The maximum number of generations which should be processed.", new IntValue(1000)));


184  Parameters.Add(new FixedValueParameter<BoolValue>(DominateOnEqualQualitiesName, "Flag which determines wether Individuals with equal quality values should be treated as dominated.", new BoolValue(false)));


185  }


186 


187  // Persistence uses this ctor to improve deserialization efficiency. If we would use the


188  // default ctor instead this would completely initialize the object (e.g. creating


189  // parameters) even though the data is later overwritten by the stored data.


190  [StorableConstructor]


191  public NSGA3(StorableConstructorFlag _) : base(_) { }


192 


193  // Each clonable item must have a cloning ctor (deep cloning, the cloner is used to handle


194  // cyclic object references). Don't forget to call the cloning ctor of the base class


195  public NSGA3(NSGA3 original, Cloner cloner) : base(original, cloner)


196  {


197  // todo: don't forget to clone storable fields


198  random = cloner.Clone(original.random);


199  solutions = new List<Solution>(original.solutions?.Select(cloner.Clone));


200  }


201 


202  public override IDeepCloneable Clone(Cloner cloner)


203  {


204  return new NSGA3(this, cloner);


205  }


206 


207  #endregion Constructors


208 


209  #region Initialization


210 


211  protected override void Initialize(CancellationToken cancellationToken)


212  {


213  base.Initialize(cancellationToken);


214 


215  PopulationSize.Value = ReferencePoint.GetNumberOfGeneratedReferencePoints(Problem.Maximization.Length);


216  InitResults();


217  InitReferencePoints();


218  InitFields();


219  Analyze();


220  }


221 


222  private void InitReferencePoints()


223  {


224  // Generate reference points and add them to results


225  ReferencePoints = ReferencePoint.GenerateReferencePoints(random, NumberOfObjectives);


226  ResultsGeneratedReferencePoints = Utility.ConvertToDoubleMatrix(ReferencePoints);


227  }


228 


229  private void InitFields()


230  {


231  random = new MersenneTwister();


232  InitSolutions();


233  }


234 


235  private void InitSolutions()


236  {


237  int minBound = 0;


238  int maxBound = 1;


239 


240  // Initialise solutions


241  solutions = new List<Solution>(PopulationSize.Value);


242  for (int i = 0; i < PopulationSize.Value; i++)


243  {


244  RealVector randomRealVector = new RealVector(Problem.Encoding.Length, random, minBound, maxBound);


245 


246  solutions.Add(new Solution(randomRealVector));


247  solutions[i].Fitness = Evaluate(solutions[i].Chromosome);


248  }


249  }


250 


251  private void InitResults()


252  {


253  Results.Add(new Result(GeneratedReferencePointsResultName, "The initially generated reference points", new DoubleMatrix()));


254  Results.Add(new Result(CurrentFrontResultName, "The Pareto Front", new DoubleMatrix()));


255  Results.Add(new Result(CurrentGenerationResultName, "The current generation", new IntValue(0)));


256  }


257 


258  #endregion Initialization


259 


260  #region Overriden Methods


261 


262  protected override void Run(CancellationToken cancellationToken)


263  {


264  while (ResultsCurrentGeneration.Value < MaximumGenerations.Value)


265  {


266  // create copies of generated reference points (to preserve the original ones for


267  // the next generation) maybe todo: use cloner?


268  ToNextGeneration(CreateCopyOfReferencePoints());


269  ResultsCurrentGeneration.Value++;


270  }


271  }


272 


273  #endregion Overriden Methods


274 


275  #region Private Methods


276 


277  private List<ReferencePoint> CreateCopyOfReferencePoints()


278  {


279  if (ReferencePoints == null) return null;


280 


281  List<ReferencePoint> referencePoints = new List<ReferencePoint>();


282  foreach (var referencePoint in ReferencePoints)


283  referencePoints.Add(new ReferencePoint(referencePoint));


284 


285  return referencePoints;


286  }


287 


288  private void Analyze()


289  {


290  ResultsSolutions = solutions.Select(s => s.Chromosome.ToArray()).ToMatrix();


291  Problem.Analyze(


292  solutions.Select(s => (Individual)new SingleEncodingIndividual(Problem.Encoding, new Scope { Variables = { new Variable(Problem.Encoding.Name, s.Chromosome) } })).ToArray(),


293  solutions.Select(s => s.Fitness).ToArray(),


294  Results,


295  random


296  );


297  }


298 


299  /// <summary>


300  /// Returns the fitness of the given <paramref name="chromosome" /> by applying the Evaluate


301  /// method of the Problem.


302  /// </summary>


303  /// <param name="chromosome"></param>


304  /// <returns></returns>


305  private double[] Evaluate(RealVector chromosome)


306  {


307  return Problem.Evaluate(new SingleEncodingIndividual(Problem.Encoding, new Scope { Variables = { new Variable(Problem.Encoding.Name, chromosome) } }), random);


308  }


309 


310  private void ToNextGeneration(List<ReferencePoint> referencePoints)


311  {


312  List<Solution> st = new List<Solution>();


313  List<Solution> qt = Mutate(Recombine(solutions));


314  List<Solution> rt = Utility.Concat(solutions, qt);


315  List<Solution> nextGeneration;


316 


317  // Do nondominated sort


318  var qualities = Utility.ToFitnessMatrix(rt);


319  // compute the pareto fronts using the DominationCalculator and discard the qualities


320  // part in the inner tuples


321  Fronts = DominationCalculator<Solution>.CalculateAllParetoFronts(rt.ToArray(), qualities, Problem.Maximization, out int[] rank, true)


322  .Select(list => new List<Solution>(list.Select(pair => pair.Item1))).ToList();


323 


324  int i = 0;


325  List<Solution> lf = null; // last front to be included


326  while (i < Fronts.Count && st.Count < PopulationSize.Value)


327  {


328  lf = Fronts[i];


329  st = Utility.Concat(st, lf);


330  i++;


331  }


332 


333  if (st.Count == PopulationSize.Value) // no selection needs to be done


334  nextGeneration = st;


335  else


336  {


337  int l = i  1;


338  nextGeneration = new List<Solution>();


339  for (int f = 0; f < l; f++)


340  nextGeneration = Utility.Concat(nextGeneration, Fronts[f]);


341  int k = PopulationSize.Value  nextGeneration.Count;


342  Normalize(st);


343  Associate(referencePoints);


344  List<Solution> solutionsToAdd = Niching(k, referencePoints);


345  nextGeneration = Utility.Concat(nextGeneration, solutionsToAdd);


346  }


347  }


348 


349  private void Normalize(List<Solution> population)


350  {


351  // Find the ideal point


352  double[] idealPoint = new double[NumberOfObjectives];


353  for (int j = 0; j < NumberOfObjectives; j++)


354  {


355  // Compute ideal point


356  idealPoint[j] = Utility.Min(s => s.Fitness[j], population);


357 


358  // Translate objectives


359  foreach (var solution in population)


360  solution.Fitness[j] = idealPoint[j];


361  }


362 


363  // Find the extreme points


364  Solution[] extremePoints = new Solution[NumberOfObjectives];


365  for (int j = 0; j < NumberOfObjectives; j++)


366  {


367  // Compute extreme points


368  double[] weights = new double[NumberOfObjectives];


369  for (int i = 0; i < NumberOfObjectives; i++) weights[i] = EPSILON;


370  weights[j] = 1;


371  double func(Solution s) => ASF(s.Fitness, weights);


372  extremePoints[j] = Utility.ArgMin(func, population);


373  }


374 


375  // Compute intercepts


376  List<double> intercepts = GetIntercepts(extremePoints.ToList());


377 


378  // Normalize objectives


379  NormalizeObjectives(intercepts, idealPoint);


380  }


381 


382  private void NormalizeObjectives(List<double> intercepts, double[] idealPoint)


383  {


384  for (int f = 0; f < Fronts.Count; f++)


385  {


386  foreach (var solution in Fronts[f])


387  {


388  for (int i = 0; i < NumberOfObjectives; i++)


389  {


390  if (Math.Abs(intercepts[i]  idealPoint[i]) > EPSILON)


391  {


392  solution.Fitness[i] = solution.Fitness[i] / (intercepts[i]  idealPoint[i]);


393  }


394  else


395  {


396  solution.Fitness[i] = solution.Fitness[i] / EPSILON;


397  }


398  }


399  }


400  }


401  }


402 


403  private void Associate(List<ReferencePoint> referencePoints)


404  {


405  for (int f = 0; f < Fronts.Count; f++)


406  {


407  foreach (var solution in Fronts[f])


408  {


409  // find reference point for which the perpendicular distance to the current


410  // solution is the lowest


411  var rpAndDist = Utility.MinArgMin(rp => GetPerpendicularDistance(rp.Values, solution.Fitness), referencePoints);


412  // associated reference point


413  var arp = rpAndDist.Item1;


414  // distance to that reference point


415  var dist = rpAndDist.Item2;


416 


417  if (f + 1 != Fronts.Count)


418  {


419  // Todo: Add member for reference point on index min_rp


420  arp.NumberOfAssociatedSolutions++;


421  }


422  else


423  {


424  // Todo: Add potential member for reference point on index min_rp


425  arp.AddPotentialAssociatedSolution(solution, dist);


426  }


427  }


428  }


429  }


430 


431  private List<Solution> Niching(int k, List<ReferencePoint> referencePoints)


432  {


433  List<Solution> solutions = new List<Solution>();


434  while (solutions.Count != k)


435  {


436  ReferencePoint min_rp = FindNicheReferencePoint(referencePoints);


437 


438  Solution chosen = SelectClusterMember(min_rp);


439  if (chosen == null)


440  {


441  referencePoints.Remove(min_rp);


442  }


443  else


444  {


445  min_rp.NumberOfAssociatedSolutions++;


446  min_rp.RemovePotentialAssociatedSolution(chosen);


447  solutions.Add(chosen);


448  }


449  }


450 


451  return solutions;


452  }


453 


454  private ReferencePoint FindNicheReferencePoint(List<ReferencePoint> referencePoints)


455  {


456  // the minimum number of associated solutions for a reference point over all reference points


457  int minNumber = Utility.Min(rp => rp.NumberOfAssociatedSolutions, referencePoints);


458 


459  // the reference points that share the number of associated solutions where that number


460  // is equal to minNumber


461  List<ReferencePoint> minAssociatedReferencePoints = new List<ReferencePoint>();


462  foreach (var referencePoint in referencePoints)


463  if (referencePoint.NumberOfAssociatedSolutions == minNumber)


464  minAssociatedReferencePoints.Add(referencePoint);


465 


466  if (minAssociatedReferencePoints.Count > 1)


467  return minAssociatedReferencePoints[random.Next(minAssociatedReferencePoints.Count)];


468  else


469  return minAssociatedReferencePoints.Single();


470  }


471 


472  private Solution SelectClusterMember(ReferencePoint referencePoint)


473  {


474  Solution chosen = null;


475  if (referencePoint.HasPotentialMember())


476  {


477  if (referencePoint.NumberOfAssociatedSolutions == 0)


478  chosen = referencePoint.FindClosestMember();


479  else


480  chosen = referencePoint.RandomMember();


481  }


482  return chosen;


483  }


484 


485  private double GetPerpendicularDistance(double[] values, double[] fitness)


486  {


487  double numerator = 0;


488  double denominator = 0;


489  for (int i = 0; i < values.Length; i++)


490  {


491  numerator += values[i] * fitness[i];


492  denominator += Math.Pow(values[i], 2);


493  }


494  double k = numerator / denominator;


495 


496  double d = 0;


497  for (int i = 0; i < values.Length; i++)


498  {


499  d += Math.Pow(k * values[i]  fitness[i], 2);


500  }


501  return Math.Sqrt(d);


502  }


503 


504  private double ASF(double[] x, double[] weight)


505  {


506  List<int> dimensions = new List<int>();


507  for (int i = 0; i < NumberOfObjectives; i++) dimensions.Add(i);


508  double f(int dim) => x[dim] / weight[dim];


509  return Utility.Max(f, dimensions);


510  }


511 


512  private List<double> GetIntercepts(List<Solution> extremePoints)


513  {


514  // Check whether there are duplicate extreme points. This might happen but the original


515  // paper does not mention how to deal with it.


516  bool duplicate = false;


517  for (int i = 0; !duplicate && i < extremePoints.Count; i++)


518  {


519  for (int j = i + 1; !duplicate && j < extremePoints.Count; j++)


520  {


521  // maybe todo: override Equals method of solution?


522  duplicate = extremePoints[i].Equals(extremePoints[j]);


523  }


524  }


525 


526  List<double> intercepts = new List<double>();


527 


528  if (duplicate)


529  { // cannot construct the unique hyperplane (this is a casual method to deal with the condition)


530  for (int f = 0; f < NumberOfObjectives; f++)


531  {


532  // extreme_points[f] stands for the individual with the largest value of


533  // objective f


534  intercepts.Add(extremePoints[f].Fitness[f]);


535  }


536  }


537  else


538  {


539  // Find the equation of the hyperplane


540  List<double> b = new List<double>(); //(pop[0].objs().size(), 1.0);


541  for (int i = 0; i < NumberOfObjectives; i++)


542  {


543  b.Add(1.0);


544  }


545 


546  List<List<double>> a = new List<List<double>>();


547  foreach (Solution s in extremePoints)


548  {


549  List<double> aux = new List<double>();


550  for (int i = 0; i < NumberOfObjectives; i++)


551  aux.Add(s.Fitness[i]);


552  a.Add(aux);


553  }


554  List<double> x = GaussianElimination(a, b);


555 


556  // Find intercepts


557  for (int f = 0; f < NumberOfObjectives; f++)


558  {


559  intercepts.Add(1.0 / x[f]);


560  }


561  }


562 


563  return intercepts;


564  }


565 


566  private List<double> GaussianElimination(List<List<double>> a, List<double> b)


567  {


568  List<double> x = new List<double>();


569 


570  int n = a.Count;


571  for (int i = 0; i < n; i++)


572  a[i].Add(b[i]);


573 


574  for (int @base = 0; @base < n  1; @base++)


575  for (int target = @base + 1; target < n; target++)


576  {


577  double ratio = a[target][@base] / a[@base][@base];


578  for (int term = 0; term < a[@base].Count; term++)


579  a[target][term] = a[target][term]  a[@base][term] * ratio;


580  }


581 


582  for (int i = 0; i < n; i++)


583  x.Add(0.0);


584 


585  for (int i = n  1; i >= 0; i)


586  {


587  for (int known = i + 1; known < n; known++)


588  a[i][n] = a[i][n]  a[i][known] * x[known];


589  x[i] = a[i][n] / a[i][i];


590  }


591 


592  return x;


593  }


594 


595  private List<Solution> Recombine(List<Solution> solutions)


596  {


597  List<Solution> childSolutions = new List<Solution>();


598 


599  for (int i = 0; i < solutions.Count; i += 2)


600  {


601  int parentIndex1 = random.Next(solutions.Count);


602  int parentIndex2 = random.Next(solutions.Count);


603  // ensure that the parents are not the same object


604  if (parentIndex1 == parentIndex2) parentIndex2 = (parentIndex2 + 1) % solutions.Count;


605  var parent1 = solutions[parentIndex1];


606  var parent2 = solutions[parentIndex2];


607 


608  // Do crossover with crossoverProbabilty == 1 in order to guarantee that a crossover happens


609  var children = SimulatedBinaryCrossover.Apply(random, Problem.Encoding.Bounds, parent1.Chromosome, parent2.Chromosome, 1);


610  Debug.Assert(children != null);


611 


612  var child1 = new Solution(children.Item1);


613  var child2 = new Solution(children.Item2);


614  child1.Fitness = Evaluate(child1.Chromosome);


615  child2.Fitness = Evaluate(child1.Chromosome);


616 


617  childSolutions.Add(child1);


618  childSolutions.Add(child2);


619  }


620 


621  return childSolutions;


622  }


623 


624  private List<Solution> Mutate(List<Solution> solutions)


625  {


626  foreach (var solution in solutions)


627  {


628  UniformOnePositionManipulator.Apply(random, solution.Chromosome, Problem.Encoding.Bounds);


629  }


630  return solutions;


631  }


632 


633  #endregion Private Methods


634  }


635  } 
