Free cookie consent management tool by TermsFeed Policy Generator

source: branches/2825-NSGA3/HeuristicLab.Algorithms.NSGA3/3.3/NSGA3.cs @ 17663

Last change on this file since 17663 was 17663, checked in by dleko, 4 years ago

#2825 Initialize population size with the correct number (according to NSGA-III paper).

File size: 25.5 KB
Line 
1using System;
2using System.Collections.Generic;
3using System.Diagnostics;
4using System.Linq;
5using System.Threading;
6using HEAL.Attic;
7using HeuristicLab.Common;
8using HeuristicLab.Core;
9using HeuristicLab.Data;
10using HeuristicLab.Encodings.RealVectorEncoding;
11using HeuristicLab.Optimization;
12using HeuristicLab.Parameters;
13using HeuristicLab.Random;
14
15namespace HeuristicLab.Algorithms.NSGA3
16{
17    /// <summary>
18    /// The Reference Point Based Non-dominated Sorting Genetic Algorithm III was introduced in Deb
19    /// et al. 2013. An Evolutionary Many-Objective Optimization Algorithm Using Reference Point
20    /// Based Non-dominated Sorting Approach. IEEE Transactions on Evolutionary Computation, 18(4),
21    /// pp. 577-601.
22    /// </summary>
23    [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.")]
24    [Creatable(Category = CreatableAttribute.Categories.PopulationBasedAlgorithms, Priority = 136)]
25    [StorableType("07C745F7-A8A3-4F99-8B2C-F97E639F9AC3")]
26    public class NSGA3 : BasicAlgorithm
27    {
28        private const double EPSILON = 10e-6; // 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 int generation;
56
57        [Storable]
58        private List<Solution> solutions;
59
60        #endregion Storable fields
61
62        #region ParameterAndResultsNames
63
64        // Parameter Names
65
66        private const string SeedName = "Seed";
67        private const string SetSeedRandomlyName = "SetSeedRandomly";
68        private const string PopulationSizeName = "PopulationSize";
69        private const string CrossoverProbabilityName = "CrossoverProbability";
70        private const string CrossoverContiguityName = "CrossoverContiguity";
71        private const string MutationProbabilityName = "MutationProbability";
72        private const string MaximumGenerationsName = "MaximumGenerations";
73        private const string DominateOnEqualQualitiesName = "DominateOnEqualQualities";
74
75        // Results Names
76
77        private const string GeneratedReferencePointsResultName = "Generated Reference Points";
78        private const string CurrentFrontResultName = "Pareto Front"; // Do not touch this
79
80        #endregion ParameterAndResultsNames
81
82        #region ParameterProperties
83
84        private IFixedValueParameter<IntValue> SeedParameter
85        {
86            get { return (IFixedValueParameter<IntValue>)Parameters[SeedName]; }
87        }
88
89        private IFixedValueParameter<BoolValue> SetSeedRandomlyParameter
90        {
91            get { return (IFixedValueParameter<BoolValue>)Parameters[SetSeedRandomlyName]; }
92        }
93
94        private IFixedValueParameter<IntValue> PopulationSizeParameter
95        {
96            get { return (IFixedValueParameter<IntValue>)Parameters[PopulationSizeName]; }
97        }
98
99        private IFixedValueParameter<PercentValue> CrossoverProbabilityParameter
100        {
101            get { return (IFixedValueParameter<PercentValue>)Parameters[CrossoverProbabilityName]; }
102        }
103
104        private IFixedValueParameter<DoubleValue> CrossoverContiguityParameter
105        {
106            get { return (IFixedValueParameter<DoubleValue>)Parameters[CrossoverContiguityName]; }
107        }
108
109        private IFixedValueParameter<PercentValue> MutationProbabilityParameter
110        {
111            get { return (IFixedValueParameter<PercentValue>)Parameters[MutationProbabilityName]; }
112        }
113
114        private IFixedValueParameter<IntValue> MaximumGenerationsParameter
115        {
116            get { return (IFixedValueParameter<IntValue>)Parameters[MaximumGenerationsName]; }
117        }
118
119        private IFixedValueParameter<BoolValue> DominateOnEqualQualitiesParameter
120        {
121            get { return (IFixedValueParameter<BoolValue>)Parameters[DominateOnEqualQualitiesName]; }
122        }
123
124        #endregion ParameterProperties
125
126        #region Properties
127
128        public IntValue Seed => SeedParameter.Value;
129
130        public BoolValue SetSeedRandomly => SetSeedRandomlyParameter.Value;
131
132        public IntValue PopulationSize => PopulationSizeParameter.Value;
133
134        public PercentValue CrossoverProbability => CrossoverProbabilityParameter.Value;
135
136        public DoubleValue CrossoverContiguity => CrossoverContiguityParameter.Value;
137
138        public PercentValue MutationProbability => MutationProbabilityParameter.Value;
139
140        public IntValue MaximumGenerations => MaximumGenerationsParameter.Value;
141
142        public BoolValue DominateOnEqualQualities => DominateOnEqualQualitiesParameter.Value;
143
144        public List<List<Solution>> Fronts { get; private set; }
145
146        public List<ReferencePoint> ReferencePoints { get; private set; }
147
148        // todo: create one property for the Generated Reference Points and one for the current
149        // generations reference points
150
151        #endregion Properties
152
153        #region ResultsProperties
154
155        public DoubleMatrix ResultsGeneratedReferencePoints
156        {
157            get { return (DoubleMatrix)Results[GeneratedReferencePointsResultName].Value; }
158            set { Results[GeneratedReferencePointsResultName].Value = value; }
159        }
160
161        public DoubleMatrix ResultsSolutions
162        {
163            get { return (DoubleMatrix)Results[CurrentFrontResultName].Value; }
164            set { Results[CurrentFrontResultName].Value = value; }
165        }
166
167        #endregion ResultsProperties
168
169        #region Constructors
170
171        public NSGA3() : base()
172        {
173            Parameters.Add(new FixedValueParameter<IntValue>(SeedName, "The random seed used to initialize the new pseudo random number generator.", new IntValue(0)));
174            Parameters.Add(new FixedValueParameter<BoolValue>(SetSeedRandomlyName, "True if the random seed should be set to a random value, otherwise false.", new BoolValue(true)));
175            Parameters.Add(new FixedValueParameter<IntValue>(PopulationSizeName, "The size of the population of Individuals.", new IntValue(ReferencePoint.GetNumberOfGeneratedReferencePoints(NumberOfObjectives))));
176            Parameters.Add(new FixedValueParameter<PercentValue>(CrossoverProbabilityName, "The probability that the crossover operator is applied on two parents.", new PercentValue(0.9)));
177            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]."));
178            Parameters.Add(new FixedValueParameter<PercentValue>(MutationProbabilityName, "The probability that the mutation operator is applied on a Individual.", new PercentValue(0.05)));
179            Parameters.Add(new FixedValueParameter<IntValue>(MaximumGenerationsName, "The maximum number of generations which should be processed.", new IntValue(1000)));
180            Parameters.Add(new FixedValueParameter<BoolValue>(DominateOnEqualQualitiesName, "Flag which determines wether Individuals with equal quality values should be treated as dominated.", new BoolValue(false)));
181        }
182
183        // Persistence uses this ctor to improve deserialization efficiency. If we would use the
184        // default ctor instead this would completely initialize the object (e.g. creating
185        // parameters) even though the data is later overwritten by the stored data.
186        [StorableConstructor]
187        public NSGA3(StorableConstructorFlag _) : base(_) { }
188
189        // Each clonable item must have a cloning ctor (deep cloning, the cloner is used to handle
190        // cyclic object references). Don't forget to call the cloning ctor of the base class
191        public NSGA3(NSGA3 original, Cloner cloner) : base(original, cloner)
192        {
193            // todo: don't forget to clone storable fields
194            random = cloner.Clone(original.random);
195            solutions = new List<Solution>(original.solutions?.Select(cloner.Clone));
196        }
197
198        public override IDeepCloneable Clone(Cloner cloner)
199        {
200            return new NSGA3(this, cloner);
201        }
202
203        #endregion Constructors
204
205        #region Initialization
206
207        protected override void Initialize(CancellationToken cancellationToken)
208        {
209            base.Initialize(cancellationToken);
210
211            InitResults();
212            InitReferencePoints();
213            InitFields();
214            Analyze();
215        }
216
217        private void InitReferencePoints()
218        {
219            // Generate reference points and add them to results
220            ReferencePoints = ReferencePoint.GenerateReferencePoints(random, NumberOfObjectives);
221            ResultsGeneratedReferencePoints = Utility.ConvertToDoubleMatrix(ReferencePoints);
222        }
223
224        private void InitFields()
225        {
226            random = new MersenneTwister();
227            generation = 0;
228            InitSolutions();
229        }
230
231        private void InitSolutions()
232        {
233            int minBound = 0;
234            int maxBound = 1;
235
236            // Initialise solutions
237            solutions = new List<Solution>(PopulationSize.Value);
238            for (int i = 0; i < PopulationSize.Value; i++)
239            {
240                RealVector randomRealVector = new RealVector(Problem.Encoding.Length, random, minBound, maxBound);
241
242                solutions.Add(new Solution(randomRealVector));
243                solutions[i].Fitness = Evaluate(solutions[i].Chromosome);
244            }
245        }
246
247        private void InitResults()
248        {
249            Results.Add(new Result(GeneratedReferencePointsResultName, "The initially generated reference points", new DoubleMatrix()));
250            Results.Add(new Result(CurrentFrontResultName, "The Pareto Front", new DoubleMatrix()));
251        }
252
253        #endregion Initialization
254
255        #region Overriden Methods
256
257        protected override void Run(CancellationToken cancellationToken)
258        {
259            while (generation != MaximumGenerations.Value)
260            {
261                // create copies of generated reference points (to preserve the original ones for
262                // the next generation) maybe todo: use cloner?
263                ToNextGeneration(CreateCopyOfReferencePoints());
264                generation++;
265            }
266        }
267
268        #endregion Overriden Methods
269
270        #region Private Methods
271
272        private List<ReferencePoint> CreateCopyOfReferencePoints()
273        {
274            if (ReferencePoints == null) return null;
275
276            List<ReferencePoint> referencePoints = new List<ReferencePoint>();
277            foreach (var referencePoint in ReferencePoints)
278                referencePoints.Add(new ReferencePoint(referencePoint));
279
280            return referencePoints;
281        }
282
283        private void Analyze()
284        {
285            ResultsSolutions = solutions.Select(s => s.Chromosome.ToArray()).ToMatrix();
286            Problem.Analyze(
287                solutions.Select(s => (Individual)new SingleEncodingIndividual(Problem.Encoding, new Scope { Variables = { new Variable(Problem.Encoding.Name, s.Chromosome) } })).ToArray(),
288                solutions.Select(s => s.Fitness).ToArray(),
289                Results,
290                random
291                );
292        }
293
294        /// <summary>
295        /// Returns the fitness of the given <paramref name="chromosome" /> by applying the Evaluate
296        /// method of the Problem.
297        /// </summary>
298        /// <param name="chromosome"></param>
299        /// <returns></returns>
300        private double[] Evaluate(RealVector chromosome)
301        {
302            return Problem.Evaluate(new SingleEncodingIndividual(Problem.Encoding, new Scope { Variables = { new Variable(Problem.Encoding.Name, chromosome) } }), random);
303        }
304
305        private void ToNextGeneration(List<ReferencePoint> referencePoints)
306        {
307            List<Solution> st = new List<Solution>();
308            List<Solution> qt = Mutate(Recombine(solutions));
309            List<Solution> rt = Utility.Concat(solutions, qt);
310            List<Solution> nextGeneration;
311
312            // Do non-dominated sort
313            var qualities = Utility.ToFitnessMatrix(rt);
314            // compute the pareto fronts using the DominationCalculator and discard the qualities
315            // part in the inner tuples
316            Fronts = DominationCalculator<Solution>.CalculateAllParetoFronts(rt.ToArray(), qualities, Problem.Maximization, out int[] rank, true)
317                .Select(list => new List<Solution>(list.Select(pair => pair.Item1))).ToList();
318
319            int i = 0;
320            List<Solution> lf = null; // last front to be included
321            while (i < Fronts.Count && st.Count < PopulationSize.Value)
322            {
323                lf = Fronts[i];
324                st = Utility.Concat(st, lf);
325                i++;
326            }
327
328            if (st.Count == PopulationSize.Value) // no selection needs to be done
329                nextGeneration = st;
330            else
331            {
332                int l = i - 1;
333                nextGeneration = new List<Solution>();
334                for (int f = 0; f < l; f++)
335                    nextGeneration = Utility.Concat(nextGeneration, Fronts[f]);
336                int k = PopulationSize.Value - nextGeneration.Count;
337                Normalize(st);
338                Associate(referencePoints);
339                List<Solution> solutionsToAdd = Niching(k, referencePoints);
340                nextGeneration = Utility.Concat(nextGeneration, solutionsToAdd);
341            }
342        }
343
344        private void Normalize(List<Solution> population)
345        {
346            // Find the ideal point
347            double[] idealPoint = new double[NumberOfObjectives];
348            for (int j = 0; j < NumberOfObjectives; j++)
349            {
350                // Compute ideal point
351                idealPoint[j] = Utility.Min(s => s.Fitness[j], population);
352
353                // Translate objectives
354                foreach (var solution in population)
355                    solution.Fitness[j] -= idealPoint[j];
356            }
357
358            // Find the extreme points
359            Solution[] extremePoints = new Solution[NumberOfObjectives];
360            for (int j = 0; j < NumberOfObjectives; j++)
361            {
362                // Compute extreme points
363                double[] weights = new double[NumberOfObjectives];
364                for (int i = 0; i < NumberOfObjectives; i++) weights[i] = EPSILON;
365                weights[j] = 1;
366                double func(Solution s) => ASF(s.Fitness, weights);
367                extremePoints[j] = Utility.ArgMin(func, population);
368            }
369
370            // Compute intercepts
371            List<double> intercepts = GetIntercepts(extremePoints.ToList());
372
373            // Normalize objectives
374            NormalizeObjectives(intercepts, idealPoint);
375        }
376
377        private void NormalizeObjectives(List<double> intercepts, double[] idealPoint)
378        {
379            for (int f = 0; f < Fronts.Count; f++)
380            {
381                foreach (var solution in Fronts[f])
382                {
383                    for (int i = 0; i < NumberOfObjectives; i++)
384                    {
385                        if (Math.Abs(intercepts[i] - idealPoint[i]) > EPSILON)
386                        {
387                            solution.Fitness[i] = solution.Fitness[i] / (intercepts[i] - idealPoint[i]);
388                        }
389                        else
390                        {
391                            solution.Fitness[i] = solution.Fitness[i] / EPSILON;
392                        }
393                    }
394                }
395            }
396        }
397
398        private void Associate(List<ReferencePoint> referencePoints)
399        {
400            for (int f = 0; f < Fronts.Count; f++)
401            {
402                foreach (var solution in Fronts[f])
403                {
404                    // find reference point for which the perpendicular distance to the current
405                    // solution is the lowest
406                    var rpAndDist = Utility.MinArgMin(rp => GetPerpendicularDistance(rp.Values, solution.Fitness), referencePoints);
407                    // associated reference point
408                    var arp = rpAndDist.Item1;
409                    // distance to that reference point
410                    var dist = rpAndDist.Item2;
411
412                    if (f + 1 != Fronts.Count)
413                    {
414                        // Todo: Add member for reference point on index min_rp
415                        arp.NumberOfAssociatedSolutions++;
416                    }
417                    else
418                    {
419                        // Todo: Add potential member for reference point on index min_rp
420                        arp.AddPotentialAssociatedSolution(solution, dist);
421                    }
422                }
423            }
424        }
425
426        private List<Solution> Niching(int k, List<ReferencePoint> referencePoints)
427        {
428            List<Solution> solutions = new List<Solution>();
429            while (solutions.Count != k)
430            {
431                ReferencePoint min_rp = FindNicheReferencePoint(referencePoints);
432
433                Solution chosen = SelectClusterMember(min_rp);
434                if (chosen == null)
435                {
436                    referencePoints.Remove(min_rp);
437                }
438                else
439                {
440                    min_rp.NumberOfAssociatedSolutions++;
441                    min_rp.RemovePotentialAssociatedSolution(chosen);
442                    solutions.Add(chosen);
443                }
444            }
445
446            return solutions;
447        }
448
449        private ReferencePoint FindNicheReferencePoint(List<ReferencePoint> referencePoints)
450        {
451            // the minimum number of associated solutions for a reference point over all reference points
452            int minNumber = Utility.Min(rp => rp.NumberOfAssociatedSolutions, referencePoints);
453
454            // the reference points that share the number of associated solutions where that number
455            // is equal to minNumber
456            List<ReferencePoint> minAssociatedReferencePoints = new List<ReferencePoint>();
457            foreach (var referencePoint in referencePoints)
458                if (referencePoint.NumberOfAssociatedSolutions == minNumber)
459                    minAssociatedReferencePoints.Add(referencePoint);
460
461            if (minAssociatedReferencePoints.Count > 1)
462                return minAssociatedReferencePoints[random.Next(minAssociatedReferencePoints.Count)];
463            else
464                return minAssociatedReferencePoints.Single();
465        }
466
467        private Solution SelectClusterMember(ReferencePoint referencePoint)
468        {
469            Solution chosen = null;
470            if (referencePoint.HasPotentialMember())
471            {
472                if (referencePoint.NumberOfAssociatedSolutions == 0)
473                    chosen = referencePoint.FindClosestMember();
474                else
475                    chosen = referencePoint.RandomMember();
476            }
477            return chosen;
478        }
479
480        private double GetPerpendicularDistance(double[] values, double[] fitness)
481        {
482            double numerator = 0;
483            double denominator = 0;
484            for (int i = 0; i < values.Length; i++)
485            {
486                numerator += values[i] * fitness[i];
487                denominator += Math.Pow(values[i], 2);
488            }
489            double k = numerator / denominator;
490
491            double d = 0;
492            for (int i = 0; i < values.Length; i++)
493            {
494                d += Math.Pow(k * values[i] - fitness[i], 2);
495            }
496            return Math.Sqrt(d);
497        }
498
499        private double ASF(double[] x, double[] weight)
500        {
501            List<int> dimensions = new List<int>();
502            for (int i = 0; i < NumberOfObjectives; i++) dimensions.Add(i);
503            double f(int dim) => x[dim] / weight[dim];
504            return Utility.Max(f, dimensions);
505        }
506
507        private List<double> GetIntercepts(List<Solution> extremePoints)
508        {
509            // Check whether there are duplicate extreme points. This might happen but the original
510            // paper does not mention how to deal with it.
511            bool duplicate = false;
512            for (int i = 0; !duplicate && i < extremePoints.Count; i++)
513            {
514                for (int j = i + 1; !duplicate && j < extremePoints.Count; j++)
515                {
516                    // maybe todo: override Equals method of solution?
517                    duplicate = extremePoints[i].Equals(extremePoints[j]);
518                }
519            }
520
521            List<double> intercepts = new List<double>();
522
523            if (duplicate)
524            { // cannot construct the unique hyperplane (this is a casual method to deal with the condition)
525                for (int f = 0; f < NumberOfObjectives; f++)
526                {
527                    // extreme_points[f] stands for the individual with the largest value of
528                    // objective f
529                    intercepts.Add(extremePoints[f].Fitness[f]);
530                }
531            }
532            else
533            {
534                // Find the equation of the hyperplane
535                List<double> b = new List<double>(); //(pop[0].objs().size(), 1.0);
536                for (int i = 0; i < NumberOfObjectives; i++)
537                {
538                    b.Add(1.0);
539                }
540
541                List<List<double>> a = new List<List<double>>();
542                foreach (Solution s in extremePoints)
543                {
544                    List<double> aux = new List<double>();
545                    for (int i = 0; i < NumberOfObjectives; i++)
546                        aux.Add(s.Fitness[i]);
547                    a.Add(aux);
548                }
549                List<double> x = GaussianElimination(a, b);
550
551                // Find intercepts
552                for (int f = 0; f < NumberOfObjectives; f++)
553                {
554                    intercepts.Add(1.0 / x[f]);
555                }
556            }
557
558            return intercepts;
559        }
560
561        private List<double> GaussianElimination(List<List<double>> a, List<double> b)
562        {
563            List<double> x = new List<double>();
564
565            int n = a.Count;
566            for (int i = 0; i < n; i++)
567                a[i].Add(b[i]);
568
569            for (int @base = 0; @base < n - 1; @base++)
570                for (int target = @base + 1; target < n; target++)
571                {
572                    double ratio = a[target][@base] / a[@base][@base];
573                    for (int term = 0; term < a[@base].Count; term++)
574                        a[target][term] = a[target][term] - a[@base][term] * ratio;
575                }
576
577            for (int i = 0; i < n; i++)
578                x.Add(0.0);
579
580            for (int i = n - 1; i >= 0; i--)
581            {
582                for (int known = i + 1; known < n; known++)
583                    a[i][n] = a[i][n] - a[i][known] * x[known];
584                x[i] = a[i][n] / a[i][i];
585            }
586
587            return x;
588        }
589
590        private List<Solution> Recombine(List<Solution> solutions)
591        {
592            List<Solution> childSolutions = new List<Solution>();
593
594            for (int i = 0; i < solutions.Count; i += 2)
595            {
596                int parentIndex1 = random.Next(solutions.Count);
597                int parentIndex2 = random.Next(solutions.Count);
598                // ensure that the parents are not the same object
599                if (parentIndex1 == parentIndex2) parentIndex2 = (parentIndex2 + 1) % solutions.Count;
600                var parent1 = solutions[parentIndex1];
601                var parent2 = solutions[parentIndex2];
602
603                // Do crossover with crossoverProbabilty == 1 in order to guarantee that a crossover happens
604                var children = SimulatedBinaryCrossover.Apply(random, Problem.Encoding.Bounds, parent1.Chromosome, parent2.Chromosome, 1);
605                Debug.Assert(children != null);
606
607                var child1 = new Solution(children.Item1);
608                var child2 = new Solution(children.Item2);
609                child1.Fitness = Evaluate(child1.Chromosome);
610                child2.Fitness = Evaluate(child1.Chromosome);
611
612                childSolutions.Add(child1);
613                childSolutions.Add(child2);
614            }
615
616            return childSolutions;
617        }
618
619        private List<Solution> Mutate(List<Solution> solutions)
620        {
621            foreach (var solution in solutions)
622            {
623                UniformOnePositionManipulator.Apply(random, solution.Chromosome, Problem.Encoding.Bounds);
624            }
625            return solutions;
626        }
627
628        #endregion Private Methods
629    }
630}
Note: See TracBrowser for help on using the repository browser.