1  using System;


2  using System.Collections.Generic;


3  using System.Linq;


4  using System.Threading;


5  using HEAL.Attic;


6  using HeuristicLab.Common;


7  using HeuristicLab.Core;


8  using HeuristicLab.Data;


9  using HeuristicLab.Encodings.RealVectorEncoding;


10  using HeuristicLab.Optimization;


11  using HeuristicLab.Parameters;


12  using HeuristicLab.Problems.TestFunctions.MultiObjective;


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  public override bool SupportsPause => true;


29 


30  #region ProblemProperties


31 


32  public override Type ProblemType


33  {


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


35  }


36 


37  public new MultiObjectiveBasicProblem<RealVectorEncoding> Problem


38  {


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


40  set { base.Problem = value; }


41  }


42 


43  public int Objectives


44  {


45  get


46  {


47  if (!(Problem is MultiObjectiveTestFunctionProblem testFunctionProblem)) throw new NotSupportedException("Only Multi Objective Test Function problems are supported");


48  return testFunctionProblem.Objectives;


49  }


50  }


51 


52  #endregion ProblemProperties


53 


54  #region Storable fields


55 


56  [Storable]


57  private IRandom random;


58 


59  [Storable]


60  private List<Solution> solutions; // maybe todo: rename to nextGeneration (see Run method)


61 


62  [Storable]


63  private List<ReferencePoint> referencePoints;


64 


65  #endregion Storable fields


66 


67  #region ParameterAndResultsNames


68 


69  // Parameter Names


70 


71  private const string PopulationSizeName = "Population Size";


72  private const string MaximumGenerationsName = "Maximum Generations";


73  private const string CrossoverProbabilityName = "Crossover Probability";


74  private const string MutationProbabilityName = "Mutation Probability";


75  private const string DominateOnEqualQualitiesName = "Dominate On Equal Qualities";


76  private const string SetSeedRandomlyName = "Set Seed Randomly";


77  private const string SeedName = "Seed";


78 


79  // Results Names


80 


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


82  private const string CurrentGenerationResultName = "Generations";


83  private const string GenerationalDistanceResultName = "Generational Distance";


84  private const string InvertedGenerationalDistanceResultName = "Inverted Generational Distance";


85  private const string HypervolumeResultName = "Hypervolume";


86  private const string BestKnownHypervolumeResultName = "Best known hypervolume";


87  private const string DifferenceToBestKnownHypervolumeResultName = "Absolute Distance to BestKnownHypervolume";


88  private const string ScatterPlotResultName = "Scatter Plot";


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


90 


91  #endregion ParameterAndResultsNames


92 


93  #region ParameterProperties


94 


95  private IFixedValueParameter<IntValue> PopulationSizeParameter


96  {


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


98  }


99 


100  private IFixedValueParameter<IntValue> MaximumGenerationsParameter


101  {


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


103  }


104 


105  private IFixedValueParameter<PercentValue> CrossoverProbabilityParameter


106  {


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


108  }


109 


110  private IFixedValueParameter<PercentValue> MutationProbabilityParameter


111  {


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


113  }


114 


115  private IFixedValueParameter<BoolValue> DominateOnEqualQualitiesParameter


116  {


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


118  }


119 


120  private IFixedValueParameter<BoolValue> SetSeedRandomlyParameter


121  {


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


123  }


124 


125  private IFixedValueParameter<IntValue> SeedParameter


126  {


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


128  }


129 


130  #endregion ParameterProperties


131 


132  #region Properties


133 


134  public IntValue PopulationSize => PopulationSizeParameter.Value;


135 


136  public IntValue MaximumGenerations => MaximumGenerationsParameter.Value;


137 


138  public PercentValue CrossoverProbability => CrossoverProbabilityParameter.Value;


139 


140  public PercentValue MutationProbability => MutationProbabilityParameter.Value;


141 


142  public BoolValue DominateOnEqualQualities => DominateOnEqualQualitiesParameter.Value;


143 


144  public BoolValue SetSeedRandomly => SetSeedRandomlyParameter.Value;


145  public IntValue Seed => SeedParameter.Value;


146 


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


148  // generations reference points


149 


150  #endregion Properties


151 


152  #region ResultsProperties


153 


154  public DoubleMatrix ResultsGeneratedReferencePoints


155  {


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


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


158  }


159 


160  public IntValue ResultsCurrentGeneration


161  {


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


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


164  }


165 


166  public DoubleValue ResultsGenerationalDistance


167  {


168  get { return (DoubleValue)Results[GenerationalDistanceResultName].Value; }


169  set { Results[GenerationalDistanceResultName].Value = value; }


170  }


171 


172  public DoubleValue ResultsInvertedGenerationalDistance


173  {


174  get { return (DoubleValue)Results[InvertedGenerationalDistanceResultName].Value; }


175  set { Results[InvertedGenerationalDistanceResultName].Value = value; }


176  }


177 


178  public DoubleValue ResultsHypervolume


179  {


180  get { return (DoubleValue)Results[HypervolumeResultName].Value; }


181  set { Results[HypervolumeResultName].Value = value; }


182  }


183 


184  public DoubleValue ResultsBestKnownHypervolume


185  {


186  get { return (DoubleValue)Results[BestKnownHypervolumeResultName].Value; }


187  set { Results[BestKnownHypervolumeResultName].Value = value; }


188  }


189 


190  public DoubleValue ResultsDifferenceToBestKnownHypervolume


191  {


192  get { return (DoubleValue)Results[DifferenceToBestKnownHypervolumeResultName].Value; }


193  set { Results[DifferenceToBestKnownHypervolumeResultName].Value = value; }


194  }


195 


196  public ParetoFrontScatterPlot ResultsScatterPlot


197  {


198  get { return (ParetoFrontScatterPlot)Results[ScatterPlotResultName].Value; }


199  set { Results[ScatterPlotResultName].Value = value; }


200  }


201 


202  public DoubleMatrix ResultsSolutions


203  {


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


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


206  }


207 


208  #endregion ResultsProperties


209 


210  #region Constructors


211 


212  public NSGA3() : base()


213  {


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


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


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


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


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


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


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


221  }


222 


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


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


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


226  [StorableConstructor]


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


228 


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


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


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


232  {


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


234  solutions = original.solutions?.Select(cloner.Clone).ToList();


235  referencePoints = original.referencePoints?.Select(cloner.Clone).ToList();


236  }


237 


238  public override IDeepCloneable Clone(Cloner cloner)


239  {


240  return new NSGA3(this, cloner);


241  }


242 


243  #endregion Constructors


244 


245  #region Initialization


246 


247  protected override void Initialize(CancellationToken cancellationToken)


248  {


249  base.Initialize(cancellationToken);


250 


251  SetParameters();


252  InitResults();


253  InitFields();


254  Analyze();


255  }


256 


257  // todo: trigger this when the problem parameters change


258  private void SetParameters()


259  {


260  // Set population size


261  int numberOfGeneratedReferencePoints = ReferencePoint.GetNumberOfGeneratedReferencePoints(Objectives);


262  if (numberOfGeneratedReferencePoints == 1) throw new NotSupportedException("The number of objectives is not supported.");


263  PopulationSize.Value = ReferencePoint.GetPopulationSizeForReferencePoints(numberOfGeneratedReferencePoints);


264 


265  // Set mutation probability


266  MutationProbability.Value = 1.0 / Problem.Encoding.Length;


267 


268  // todo: Set MaximumGenerations in Problem definition


269  }


270 


271  private void InitResults()


272  {


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


274  Results.Add(new Result(CurrentGenerationResultName, "The current generation", new IntValue(1)));


275  Results.Add(new Result(GenerationalDistanceResultName, "The generational distance to an optimal pareto front defined in the Problem", new DoubleValue(double.NaN)));


276  Results.Add(new Result(InvertedGenerationalDistanceResultName, "The inverted generational distance to an optimal pareto front defined in the Problem", new DoubleValue(double.NaN)));


277  Results.Add(new Result(HypervolumeResultName, "The hypervolume of the current front considering the Reference point defined in the Problem", new DoubleValue(double.NaN)));


278  Results.Add(new Result(BestKnownHypervolumeResultName, "The best known hypervolume considering the Reference point defined in the Problem", new DoubleValue(double.NaN)));


279  Results.Add(new Result(DifferenceToBestKnownHypervolumeResultName, "The difference between the current and the best known hypervolume", new DoubleValue(double.NaN)));


280  Results.Add(new Result(ScatterPlotResultName, "A scatterplot displaying the evaluated solutions and (if available) the analytically optimal front", new ParetoFrontScatterPlot()));


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


282 


283  if (!(Problem is MultiObjectiveTestFunctionProblem problem)) return;


284  // todo: add BestKnownFront parameter


285  ResultsScatterPlot = new ParetoFrontScatterPlot(new double[0][], new double[0][], problem.BestKnownFront.ToJaggedArray(), problem.Objectives, problem.ProblemSize);


286  if (problem.BestKnownFront == null) return;


287  ResultsBestKnownHypervolume = new DoubleValue(Hypervolume.Calculate(problem.BestKnownFront.ToJaggedArray(), problem.ReferencePoint.CloneAsArray(), problem.Maximization));


288  }


289 


290  private void InitFields()


291  {


292  random = new MersenneTwister();


293  if (!SetSeedRandomly.Value)


294  random.Reset(Seed.Value);


295 


296  solutions = GetInitialPopulation();


297  referencePoints = ReferencePoint.GenerateReferencePoints(random, Objectives);


298  ResultsGeneratedReferencePoints = Utility.ConvertToDoubleMatrix(referencePoints);


299  }


300 


301  private List<Solution> GetInitialPopulation()


302  {


303  var problem = Problem as MultiObjectiveTestFunctionProblem;


304  if (problem.Bounds.Rows != 1) throw new Exception();


305 


306  // Initialise solutions


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


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


309  {


310  double minBound = problem.Bounds[0, 0];


311  double maxBound = problem.Bounds[0, 1];


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


313  var solution = new Solution(randomRealVector);


314  solution.Fitness = Evaluate(solution.Chromosome);


315  solutions.Add(solution);


316  }


317 


318  return solutions;


319  }


320 


321  #endregion Initialization


322 


323  #region Overriden Methods


324 


325  protected override void Run(CancellationToken cancellationToken)


326  {


327  try


328  {


329  while (ResultsCurrentGeneration.Value < MaximumGenerations.Value)


330  {


331  // todo: make parameter out of this


332  List<Solution> qt = Mutate(Recombine(solutions), 30);


333  foreach (var solution in qt)


334  solution.Fitness = Evaluate(solution.Chromosome);


335 


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


337 


338  // todo: remove this check


339  for (int i = 0; i < rt.Count / 2; i++)


340  if (!solutions.Contains(rt[i])) throw new Exception($"This should never happen: !solutions.Contains(rt[{i}])");


341 


342  solutions = NSGA3Selection.SelectSolutionsForNextGeneration(rt, GetCopyOfReferencePoints(), Problem.Maximization, PopulationSize.Value, random);


343 


344  ResultsCurrentGeneration.Value++;


345  cancellationToken.ThrowIfCancellationRequested();


346  }


347  }


348  catch (OperationCanceledException ex)


349  {


350  throw new OperationCanceledException("Optimization process was cancelled.", ex);


351  }


352  catch (Exception ex)


353  {


354  throw new Exception($"Failed in generation {ResultsCurrentGeneration}.", ex);


355  }


356  finally


357  {


358  Analyze();


359  }


360  }


361 


362  #endregion Overriden Methods


363 


364  #region Private Methods


365 


366  private List<ReferencePoint> GetCopyOfReferencePoints()


367  {


368  if (referencePoints == null) return null;


369 


370  return referencePoints.Select(rp => (ReferencePoint)rp.Clone()).ToList();


371  }


372 


373  private void Analyze()


374  {


375  ResultsScatterPlot = new ParetoFrontScatterPlot(solutions.Select(x => x.Fitness).ToArray(), solutions.Select(x => x.Chromosome.ToArray()).ToArray(), ResultsScatterPlot.ParetoFront, ResultsScatterPlot.Objectives, ResultsScatterPlot.ProblemSize);


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


377 


378  if (!(Problem is MultiObjectiveTestFunctionProblem problem)) return;


379 


380  // Indicators


381  ResultsGenerationalDistance = new DoubleValue(problem.BestKnownFront != null ? GenerationalDistance.Calculate(solutions.Select(s => s.Fitness), problem.BestKnownFront.ToJaggedArray(), 1) : double.NaN);


382  ResultsInvertedGenerationalDistance = new DoubleValue(problem.BestKnownFront != null ? InvertedGenerationalDistance.Calculate(solutions.Select(s => s.Fitness), problem.BestKnownFront.ToJaggedArray(), 1) : double.NaN);


383 


384  var front = NonDominatedSelect.GetDominatingVectors(solutions.Select(x => x.Fitness), problem.ReferencePoint.CloneAsArray(), Problem.Maximization, true).ToArray();


385  if (front.Length == 0) return;


386  ResultsHypervolume = new DoubleValue(Hypervolume.Calculate(front, problem.ReferencePoint.CloneAsArray(), problem.Maximization));


387  ResultsDifferenceToBestKnownHypervolume = new DoubleValue(ResultsBestKnownHypervolume.Value  ResultsHypervolume.Value);


388  }


389 


390  /// <summary>


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


392  /// method of the Problem.


393  /// </summary>


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


395  /// <returns></returns>


396  private double[] Evaluate(RealVector chromosome)


397  {


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


399  }


400 


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


402  {


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


404 


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


406  {


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


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


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


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


411  var parent1 = solutions[parentIndex1];


412  var parent2 = solutions[parentIndex2];


413 


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


415  var children = SimulatedBinaryCrossover.Apply(random,


416  Problem.Encoding.Bounds, parent1.Chromosome, parent2.Chromosome, CrossoverProbability.Value);


417 


418  childSolutions.Add(new Solution(children.Item1));


419  childSolutions.Add(new Solution(children.Item2));


420  }


421 


422  return childSolutions;


423  }


424 


425  private List<Solution> Mutate(List<Solution> solutions, double eta)


426  {


427  foreach (var solution in solutions)


428  {


429  for (int i = 0; i < solution.Chromosome.Length; i++)


430  {


431  if (random.NextDouble() > MutationProbability.Value) continue;


432 


433  double y = solution.Chromosome[i];


434  double lb;


435  double ub;


436  var problem = Problem as MultiObjectiveTestFunctionProblem;


437  if (problem.Bounds.Rows == 1) lb = problem.Bounds[0, 0];


438  else lb = problem.Bounds[i, 0];


439  if (problem.Bounds.Rows == 1) ub = problem.Bounds[0, 1];


440  else ub = problem.Bounds[i, 1];


441 


442  double delta1 = (y  lb) / (ub  lb);


443  double delta2 = (ub  y) / (ub  lb);


444 


445  double mut_pow = 1.0 / (eta + 1.0);


446 


447  double rnd = random.NextDouble();


448  double deltaq;


449  if (rnd <= 0.5)


450  {


451  double xy = 1.0  delta1;


452  double val = 2.0 * rnd + (1.0  2.0 * rnd) * (Math.Pow(xy, (eta + 1.0)));


453  deltaq = Math.Pow(val, mut_pow)  1.0;


454  }


455  else


456  {


457  double xy = 1.0  delta2;


458  double val = 2.0 * (1.0  rnd) + 2.0 * (rnd  0.5) * (Math.Pow(xy, (eta + 1.0)));


459  deltaq = 1.0  (Math.Pow(val, mut_pow));


460  }


461 


462  y += deltaq * (ub  lb);


463  y = Math.Min(ub, Math.Max(lb, y));


464 


465  solution.Chromosome[i] = y;


466  }


467  }


468  return solutions;


469  }


470 


471  #endregion Private Methods


472  }


473  } 
