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 NumberOfObjectives


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 ScatterPlotResultName = "Scatter Plot";


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


87 


88  #endregion ParameterAndResultsNames


89 


90  #region ParameterProperties


91 


92  private IFixedValueParameter<IntValue> PopulationSizeParameter


93  {


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


95  }


96 


97  private IFixedValueParameter<IntValue> MaximumGenerationsParameter


98  {


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


100  }


101 


102  private IFixedValueParameter<PercentValue> CrossoverProbabilityParameter


103  {


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


105  }


106 


107  private IFixedValueParameter<PercentValue> MutationProbabilityParameter


108  {


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


110  }


111 


112  private IFixedValueParameter<BoolValue> DominateOnEqualQualitiesParameter


113  {


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


115  }


116 


117  private IFixedValueParameter<BoolValue> SetSeedRandomlyParameter


118  {


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


120  }


121 


122  private IFixedValueParameter<IntValue> SeedParameter


123  {


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


125  }


126 


127  #endregion ParameterProperties


128 


129  #region Properties


130 


131  public IntValue PopulationSize => PopulationSizeParameter.Value;


132 


133  public IntValue MaximumGenerations => MaximumGenerationsParameter.Value;


134 


135  public PercentValue CrossoverProbability => CrossoverProbabilityParameter.Value;


136 


137  public PercentValue MutationProbability => MutationProbabilityParameter.Value;


138 


139  public BoolValue DominateOnEqualQualities => DominateOnEqualQualitiesParameter.Value;


140 


141  public BoolValue SetSeedRandomly => SetSeedRandomlyParameter.Value;


142  public IntValue Seed => SeedParameter.Value;


143 


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


145  // generations reference points


146 


147  #endregion Properties


148 


149  #region ResultsProperties


150 


151  public DoubleMatrix ResultsGeneratedReferencePoints


152  {


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


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


155  }


156 


157  public IntValue ResultsCurrentGeneration


158  {


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


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


161  }


162 


163  public DoubleValue ResultsGenerationalDistance


164  {


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


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


167  }


168 


169  public DoubleValue ResultsInvertedGenerationalDistance


170  {


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


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


173  }


174 


175  public ParetoFrontScatterPlot ResultsScatterPlot


176  {


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


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


179  }


180 


181  public DoubleMatrix ResultsSolutions


182  {


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


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


185  }


186 


187  #endregion ResultsProperties


188 


189  #region Constructors


190 


191  public NSGA3() : base()


192  {


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


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


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


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


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


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


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


200  }


201 


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


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


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


205  [StorableConstructor]


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


207 


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


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


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


211  {


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


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


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


215  referencePoints = original.referencePoints?.Select(r =>


216  {


217  var refPoint = new ReferencePoint(random, r.NumberOfDimensions);


218  r.Values.CopyTo(refPoint.Values, 0);


219  return refPoint;


220  }).ToList();


221  }


222 


223  public override IDeepCloneable Clone(Cloner cloner)


224  {


225  return new NSGA3(this, cloner);


226  }


227 


228  #endregion Constructors


229 


230  #region Initialization


231 


232  protected override void Initialize(CancellationToken cancellationToken)


233  {


234  base.Initialize(cancellationToken);


235 


236  SetParameters();


237  InitResults();


238  InitFields();


239  Analyze();


240  }


241 


242  private void SetParameters()


243  {


244  // Set population size


245  int numberOfGeneratedReferencePoints = ReferencePoint.GetNumberOfGeneratedReferencePoints(NumberOfObjectives);


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


247 


248  // Set mutation probability


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


250  }


251 


252  private void InitResults()


253  {


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


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


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


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


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


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


260 


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


262  // todo: add BestKnownFront parameter


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


264  }


265 


266  private void InitFields()


267  {


268  random = new MersenneTwister();


269  solutions = GetInitialPopulation();


270  InitReferencePoints();


271  }


272 


273  private void InitReferencePoints()


274  {


275  // Generate reference points and add them to results


276  referencePoints = ReferencePoint.GenerateReferencePoints(random, NumberOfObjectives);


277  ResultsGeneratedReferencePoints = Utility.ConvertToDoubleMatrix(referencePoints);


278  }


279 


280  private List<Solution> GetInitialPopulation()


281  {


282  var problem = Problem as MultiObjectiveTestFunctionProblem;


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


284 


285  // Initialise solutions


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


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


288  {


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


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


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


292  var solution = new Solution(randomRealVector);


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


294  solutions.Add(solution);


295  }


296 


297  return solutions;


298  }


299 


300  #endregion Initialization


301 


302  #region Overriden Methods


303 


304  protected override void Run(CancellationToken cancellationToken)


305  {


306  while (ResultsCurrentGeneration.Value < MaximumGenerations.Value)


307  {


308  try


309  {


310  // todo: make parameter out of this


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


312  foreach (var solution in qt)


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


314 


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


316 


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


318 


319  ResultsCurrentGeneration.Value++;


320  Analyze();


321  cancellationToken.ThrowIfCancellationRequested();


322  }


323  catch (OperationCanceledException ex)


324  {


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


326  }


327  catch (Exception ex)


328  {


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


330  }


331  finally


332  {


333  Analyze();


334  }


335  }


336  }


337 


338  #endregion Overriden Methods


339 


340  #region Private Methods


341 


342  private List<ReferencePoint> GetCopyOfReferencePoints()


343  {


344  if (referencePoints == null) return null;


345 


346  List<ReferencePoint> referencePointsCopy = new List<ReferencePoint>();


347  foreach (var referencePoint in referencePoints)


348  referencePointsCopy.Add(new ReferencePoint(referencePoint));


349 


350  return referencePointsCopy;


351  }


352 


353  private void Analyze()


354  {


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


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


357 


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


359 


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


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


362 


363  Problem.Analyze(


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


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


366  Results,


367  random


368  );


369  }


370 


371  /// <summary>


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


373  /// method of the Problem.


374  /// </summary>


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


376  /// <returns></returns>


377  private double[] Evaluate(RealVector chromosome)


378  {


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


380  }


381 


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


383  {


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


385 


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


387  {


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


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


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


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


392  var parent1 = solutions[parentIndex1];


393  var parent2 = solutions[parentIndex2];


394 


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


396  var children = SimulatedBinaryCrossover.Apply(random,


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


398 


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


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


401  }


402 


403  return childSolutions;


404  }


405 


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


407  {


408  foreach (var solution in solutions)


409  {


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


411  {


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


413 


414  double y = solution.Chromosome[i];


415  double lb;


416  double ub;


417  var problem = Problem as MultiObjectiveTestFunctionProblem;


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


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


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


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


422 


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


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


425 


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


427 


428  double rnd = random.NextDouble();


429  double deltaq;


430  if (rnd <= 0.5)


431  {


432  double xy = 1.0  delta1;


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


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


435  }


436  else


437  {


438  double xy = 1.0  delta2;


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


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


441  }


442 


443  y += deltaq * (ub  lb);


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


445 


446  solution.Chromosome[i] = y;


447  }


448  }


449  return solutions;


450  }


451 


452  #endregion Private Methods


453  }


454  } 
