using System; using System.Collections.Generic; using System.Linq; using System.Threading; using HEAL.Attic; using HeuristicLab.Common; using HeuristicLab.Core; using HeuristicLab.Data; using HeuristicLab.Encodings.RealVectorEncoding; using HeuristicLab.Optimization; using HeuristicLab.Parameters; using HeuristicLab.Problems.TestFunctions.MultiObjective; using HeuristicLab.Random; namespace HeuristicLab.Algorithms.NSGA3 { /// /// The Reference Point Based Non-dominated Sorting Genetic Algorithm III was introduced in Deb /// et al. 2013. An Evolutionary Many-Objective Optimization Algorithm Using Reference Point /// Based Non-dominated Sorting Approach. IEEE Transactions on Evolutionary Computation, 18(4), /// pp. 577-601. /// [Item("NSGA-III", "The Reference Point Based Non-dominated Sorting Genetic Algorithm III was introduced in Deb et al. 2013. An Evolutionary Many-Objective Optimization Algorithm Using Reference Point Based Non-dominated Sorting Approach. IEEE Transactions on Evolutionary Computation, 18(4), pp. 577-601.")] [Creatable(Category = CreatableAttribute.Categories.PopulationBasedAlgorithms, Priority = 136)] [StorableType("07C745F7-A8A3-4F99-8B2C-F97E639F9AC3")] public class NSGA3 : BasicAlgorithm { public override bool SupportsPause => true; #region ProblemProperties public override Type ProblemType { get { return typeof(MultiObjectiveBasicProblem); } } public new MultiObjectiveBasicProblem Problem { get { return (MultiObjectiveBasicProblem)base.Problem; } set { base.Problem = value; } } public int NumberOfObjectives { get { if (!(Problem is MultiObjectiveTestFunctionProblem testFunctionProblem)) throw new NotSupportedException("Only Multi Objective Test Function problems are supported"); return testFunctionProblem.Objectives; } } #endregion ProblemProperties #region Storable fields [Storable] private IRandom random; [Storable] private List solutions; // maybe todo: rename to nextGeneration (see Run method) [Storable] private List referencePoints; #endregion Storable fields #region ParameterAndResultsNames // Parameter Names private const string PopulationSizeName = "Population Size"; private const string MaximumGenerationsName = "Maximum Generations"; private const string CrossoverProbabilityName = "Crossover Probability"; private const string MutationProbabilityName = "Mutation Probability"; private const string DominateOnEqualQualitiesName = "Dominate On Equal Qualities"; private const string SetSeedRandomlyName = "Set Seed Randomly"; private const string SeedName = "Seed"; // Results Names private const string GeneratedReferencePointsResultName = "Generated Reference Points"; private const string CurrentGenerationResultName = "Generations"; private const string GenerationalDistanceResultName = "Generational Distance"; private const string InvertedGenerationalDistanceResultName = "Inverted Generational Distance"; private const string ScatterPlotResultName = "Scatter Plot"; private const string CurrentFrontResultName = "Pareto Front"; // Do not touch this #endregion ParameterAndResultsNames #region ParameterProperties private IFixedValueParameter PopulationSizeParameter { get { return (IFixedValueParameter)Parameters[PopulationSizeName]; } } private IFixedValueParameter MaximumGenerationsParameter { get { return (IFixedValueParameter)Parameters[MaximumGenerationsName]; } } private IFixedValueParameter CrossoverProbabilityParameter { get { return (IFixedValueParameter)Parameters[CrossoverProbabilityName]; } } private IFixedValueParameter MutationProbabilityParameter { get { return (IFixedValueParameter)Parameters[MutationProbabilityName]; } } private IFixedValueParameter DominateOnEqualQualitiesParameter { get { return (IFixedValueParameter)Parameters[DominateOnEqualQualitiesName]; } } private IFixedValueParameter SetSeedRandomlyParameter { get { return (IFixedValueParameter)Parameters[SetSeedRandomlyName]; } } private IFixedValueParameter SeedParameter { get { return (IFixedValueParameter)Parameters[SeedName]; } } #endregion ParameterProperties #region Properties public IntValue PopulationSize => PopulationSizeParameter.Value; public IntValue MaximumGenerations => MaximumGenerationsParameter.Value; public PercentValue CrossoverProbability => CrossoverProbabilityParameter.Value; public PercentValue MutationProbability => MutationProbabilityParameter.Value; public BoolValue DominateOnEqualQualities => DominateOnEqualQualitiesParameter.Value; public BoolValue SetSeedRandomly => SetSeedRandomlyParameter.Value; public IntValue Seed => SeedParameter.Value; // todo: create one property for the Generated Reference Points and one for the current // generations reference points #endregion Properties #region ResultsProperties public DoubleMatrix ResultsGeneratedReferencePoints { get { return (DoubleMatrix)Results[GeneratedReferencePointsResultName].Value; } set { Results[GeneratedReferencePointsResultName].Value = value; } } public IntValue ResultsCurrentGeneration { get { return (IntValue)Results[CurrentGenerationResultName].Value; } set { Results[CurrentGenerationResultName].Value = value; } } public DoubleValue ResultsGenerationalDistance { get { return (DoubleValue)Results[GenerationalDistanceResultName].Value; } set { Results[GenerationalDistanceResultName].Value = value; } } public DoubleValue ResultsInvertedGenerationalDistance { get { return (DoubleValue)Results[InvertedGenerationalDistanceResultName].Value; } set { Results[InvertedGenerationalDistanceResultName].Value = value; } } public ParetoFrontScatterPlot ResultsScatterPlot { get { return (ParetoFrontScatterPlot)Results[ScatterPlotResultName].Value; } set { Results[ScatterPlotResultName].Value = value; } } public DoubleMatrix ResultsSolutions { get { return (DoubleMatrix)Results[CurrentFrontResultName].Value; } set { Results[CurrentFrontResultName].Value = value; } } #endregion ResultsProperties #region Constructors public NSGA3() : base() { Parameters.Add(new FixedValueParameter(PopulationSizeName, "The size of the population of Individuals.", new IntValue(200))); Parameters.Add(new FixedValueParameter(MaximumGenerationsName, "The maximum number of generations which should be processed.", new IntValue(1000))); Parameters.Add(new FixedValueParameter(CrossoverProbabilityName, "The probability that the crossover operator is applied on two parents.", new PercentValue(1.0))); Parameters.Add(new FixedValueParameter(MutationProbabilityName, "The probability that the mutation operator is applied on a Individual.", new PercentValue(0.05))); Parameters.Add(new FixedValueParameter(DominateOnEqualQualitiesName, "Flag which determines wether Individuals with equal quality values should be treated as dominated.", new BoolValue(false))); Parameters.Add(new FixedValueParameter(SetSeedRandomlyName, "True if the random seed should be set to a random value, otherwise false.", new BoolValue(true))); Parameters.Add(new FixedValueParameter(SeedName, "The random seed used to initialize the new pseudo random number generator.", new IntValue(0))); } // Persistence uses this ctor to improve deserialization efficiency. If we would use the // default ctor instead this would completely initialize the object (e.g. creating // parameters) even though the data is later overwritten by the stored data. [StorableConstructor] public NSGA3(StorableConstructorFlag _) : base(_) { } // Each clonable item must have a cloning ctor (deep cloning, the cloner is used to handle // cyclic object references). Don't forget to call the cloning ctor of the base class public NSGA3(NSGA3 original, Cloner cloner) : base(original, cloner) { // todo: don't forget to clone storable fields random = cloner.Clone(original.random); solutions = original.solutions?.Select(cloner.Clone).ToList(); referencePoints = original.referencePoints?.Select(r => { var refPoint = new ReferencePoint(random, r.NumberOfDimensions); r.Values.CopyTo(refPoint.Values, 0); return refPoint; }).ToList(); } public override IDeepCloneable Clone(Cloner cloner) { return new NSGA3(this, cloner); } #endregion Constructors #region Initialization protected override void Initialize(CancellationToken cancellationToken) { base.Initialize(cancellationToken); SetParameters(); InitResults(); InitFields(); Analyze(); } private void SetParameters() { // Set population size int numberOfGeneratedReferencePoints = ReferencePoint.GetNumberOfGeneratedReferencePoints(NumberOfObjectives); PopulationSize.Value = ReferencePoint.GetPopulationSizeForReferencePoints(numberOfGeneratedReferencePoints); // Set mutation probability MutationProbability.Value = 1.0 / Problem.Encoding.Length; } private void InitResults() { Results.Add(new Result(GeneratedReferencePointsResultName, "The initially generated reference points", new DoubleMatrix())); Results.Add(new Result(CurrentGenerationResultName, "The current generation", new IntValue(1))); Results.Add(new Result(GenerationalDistanceResultName, "The generational distance to an optimal pareto front defined in the Problem", new DoubleValue(double.NaN))); Results.Add(new Result(InvertedGenerationalDistanceResultName, "The inverted generational distance to an optimal pareto front defined in the Problem", new DoubleValue(double.NaN))); Results.Add(new Result(ScatterPlotResultName, "A scatterplot displaying the evaluated solutions and (if available) the analytically optimal front", new ParetoFrontScatterPlot())); Results.Add(new Result(CurrentFrontResultName, "The Pareto Front", new DoubleMatrix())); if (!(Problem is MultiObjectiveTestFunctionProblem problem)) return; // todo: add BestKnownFront parameter ResultsScatterPlot = new ParetoFrontScatterPlot(new double[0][], new double[0][], problem.BestKnownFront.ToJaggedArray(), problem.Objectives, problem.ProblemSize); } private void InitFields() { random = new MersenneTwister(); solutions = GetInitialPopulation(); InitReferencePoints(); } private void InitReferencePoints() { // Generate reference points and add them to results referencePoints = ReferencePoint.GenerateReferencePoints(random, NumberOfObjectives); ResultsGeneratedReferencePoints = Utility.ConvertToDoubleMatrix(referencePoints); } private List GetInitialPopulation() { var problem = Problem as MultiObjectiveTestFunctionProblem; if (problem.Bounds.Rows != 1) throw new Exception(); // Initialise solutions List solutions = new List(PopulationSize.Value); for (int i = 0; i < PopulationSize.Value; i++) { double minBound = problem.Bounds[0, 0]; double maxBound = problem.Bounds[0, 1]; RealVector randomRealVector = new RealVector(Problem.Encoding.Length, random, minBound, maxBound); var solution = new Solution(randomRealVector); solution.Fitness = Evaluate(solution.Chromosome); solutions.Add(solution); } return solutions; } #endregion Initialization #region Overriden Methods protected override void Run(CancellationToken cancellationToken) { while (ResultsCurrentGeneration.Value < MaximumGenerations.Value) { try { // todo: make parameter out of this List qt = Mutate(Recombine(solutions), 30); foreach (var solution in qt) solution.Fitness = Evaluate(solution.Chromosome); List rt = Utility.Concat(solutions, qt); solutions = NSGA3Selection.SelectSolutionsForNextGeneration(rt, GetCopyOfReferencePoints(), Problem.Maximization, PopulationSize.Value, random); ResultsCurrentGeneration.Value++; Analyze(); cancellationToken.ThrowIfCancellationRequested(); } catch (OperationCanceledException ex) { throw new OperationCanceledException("Optimization process was cancelled.", ex); } catch (Exception ex) { throw new Exception($"Failed in generation {ResultsCurrentGeneration}.", ex); } finally { Analyze(); } } } #endregion Overriden Methods #region Private Methods private List GetCopyOfReferencePoints() { if (referencePoints == null) return null; List referencePointsCopy = new List(); foreach (var referencePoint in referencePoints) referencePointsCopy.Add(new ReferencePoint(referencePoint)); return referencePointsCopy; } private void Analyze() { ResultsScatterPlot = new ParetoFrontScatterPlot(solutions.Select(x => x.Fitness).ToArray(), solutions.Select(x => x.Chromosome.ToArray()).ToArray(), ResultsScatterPlot.ParetoFront, ResultsScatterPlot.Objectives, ResultsScatterPlot.ProblemSize); ResultsSolutions = solutions.Select(s => s.Chromosome.ToArray()).ToMatrix(); if (!(Problem is MultiObjectiveTestFunctionProblem problem)) return; ResultsGenerationalDistance = new DoubleValue(problem.BestKnownFront != null ? GenerationalDistance.Calculate(solutions.Select(s => s.Fitness), problem.BestKnownFront.ToJaggedArray(), 1) : double.NaN); ResultsInvertedGenerationalDistance = new DoubleValue(problem.BestKnownFront != null ? InvertedGenerationalDistance.Calculate(solutions.Select(s => s.Fitness), problem.BestKnownFront.ToJaggedArray(), 1) : double.NaN); Problem.Analyze( solutions.Select(s => (Individual)new SingleEncodingIndividual(Problem.Encoding, new Scope { Variables = { new Variable(Problem.Encoding.Name, s.Chromosome) } })).ToArray(), solutions.Select(s => s.Fitness).ToArray(), Results, random ); } /// /// Returns the fitness of the given by applying the Evaluate /// method of the Problem. /// /// /// private double[] Evaluate(RealVector chromosome) { return Problem.Evaluate(new SingleEncodingIndividual(Problem.Encoding, new Scope { Variables = { new Variable(Problem.Encoding.Name, chromosome) } }), random); } private List Recombine(List solutions) { List childSolutions = new List(); for (int i = 0; i < solutions.Count; i += 2) { int parentIndex1 = random.Next(solutions.Count); int parentIndex2 = random.Next(solutions.Count); // ensure that the parents are not the same object if (parentIndex1 == parentIndex2) parentIndex2 = (parentIndex2 + 1) % solutions.Count; var parent1 = solutions[parentIndex1]; var parent2 = solutions[parentIndex2]; // Do crossover with crossoverProbabilty == 1 in order to guarantee that a crossover happens var children = SimulatedBinaryCrossover.Apply(random, Problem.Encoding.Bounds, parent1.Chromosome, parent2.Chromosome, CrossoverProbability.Value); childSolutions.Add(new Solution(children.Item1)); childSolutions.Add(new Solution(children.Item2)); } return childSolutions; } private List Mutate(List solutions, double eta) { foreach (var solution in solutions) { for (int i = 0; i < solution.Chromosome.Length; i++) { if (random.NextDouble() > MutationProbability.Value) continue; double y = solution.Chromosome[i]; double lb; double ub; var problem = Problem as MultiObjectiveTestFunctionProblem; if (problem.Bounds.Rows == 1) lb = problem.Bounds[0, 0]; else lb = problem.Bounds[i, 0]; if (problem.Bounds.Rows == 1) ub = problem.Bounds[0, 1]; else ub = problem.Bounds[i, 1]; double delta1 = (y - lb) / (ub - lb); double delta2 = (ub - y) / (ub - lb); double mut_pow = 1.0 / (eta + 1.0); double rnd = random.NextDouble(); double deltaq; if (rnd <= 0.5) { double xy = 1.0 - delta1; double val = 2.0 * rnd + (1.0 - 2.0 * rnd) * (Math.Pow(xy, (eta + 1.0))); deltaq = Math.Pow(val, mut_pow) - 1.0; } else { double xy = 1.0 - delta2; double val = 2.0 * (1.0 - rnd) + 2.0 * (rnd - 0.5) * (Math.Pow(xy, (eta + 1.0))); deltaq = 1.0 - (Math.Pow(val, mut_pow)); } y += deltaq * (ub - lb); y = Math.Min(ub, Math.Max(lb, y)); solution.Chromosome[i] = y; } } return solutions; } #endregion Private Methods } }