diff --git a/GaSchedule.Algorithm/Cso.cs b/GaSchedule.Algorithm/Cso.cs index 64e1ca0..12a781e 100644 --- a/GaSchedule.Algorithm/Cso.cs +++ b/GaSchedule.Algorithm/Cso.cs @@ -17,9 +17,10 @@ public class Cso : NsgaIII where T : Chromosome { private int _max_iterations = 5000; private int _chromlen; - private double _pa, _beta, _σu, _σv; - private float[] _gBestScore = null; + private double _pa; + private float[] _gBest = null; private float[][] _current_position = null; + private LévyFlights _lf; // Initializes Cuckoo Search Optimization public Cso(T prototype, int numberOfCrossoverPoints = 2, int mutationSize = 2, float crossoverProbability = 80, float mutationProbability = 3) : base(prototype, numberOfCrossoverPoints, mutationSize, crossoverProbability, mutationProbability) @@ -29,12 +30,6 @@ public Cso(T prototype, int numberOfCrossoverPoints = 2, int mutationSize = 2, f _populationSize = 5; _pa = .25; - _beta = 1.5; - - var num = Gamma(1 + _beta) * Math.Sin(Math.PI * _beta / 2); - var den = Gamma((1 + _beta) / 2) * _beta * Math.Pow(2, (_beta - 1) / 2); - _σu = Math.Pow(num / den, 1 / _beta); - _σv = 1; } static E[][] CreateArray(int rows, int cols) @@ -45,25 +40,6 @@ static E[][] CreateArray(int rows, int cols) return array; } - - static double Gamma(double z) - { - if (z < 0.5) - return Math.PI / Math.Sin(Math.PI * z) / Gamma(1.0 - z); - - // Lanczos approximation g=5, n=7 - var coef = new double[7] { 1.000000000190015, 76.18009172947146, -86.50532032941677, - 24.01409824083091, -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5 }; - - var zz = z - 1.0; - var b = zz + 5.5; // g + 0.5 - var sum = coef[0]; - for (int i = 1; i < coef.Length; ++i) - sum += coef[i] / (zz + i); - - var LogSqrtTwoPi = 0.91893853320467274178; - return Math.Exp(LogSqrtTwoPi + Math.Log(sum) - b + Math.Log(b) * (zz + 0.5)); - } protected override void Initialize(List population) { @@ -76,48 +52,13 @@ protected override void Initialize(List population) if(i < 1) { _chromlen = positions.Count; _current_position = CreateArray(_populationSize, _chromlen); + _lf = new LévyFlights(_chromlen); } } } - private float[] Optimum(float[] localVal, T chromosome) - { - var localBest = _prototype.MakeEmptyFromPrototype(); - localBest.UpdatePositions(localVal); - - if(localBest.Dominates(chromosome)) { - chromosome.UpdatePositions(localVal); - return localVal; - } - - var positions = new float[_chromlen]; - chromosome.ExtractPositions(positions); - return positions; - } - - private void UpdatePosition1(List population) - { - var current_position = _current_position.ToArray(); - for(int i = 0; i < _populationSize; ++i) { - double u = Configuration.NextGaussian() * _σu; - double v = Configuration.NextGaussian() * _σv; - double S = u / Math.Pow(Math.Abs(v), 1 / _beta); - - if(_gBestScore == null) { - _gBestScore = new float[_chromlen]; - population[i].ExtractPositions(_gBestScore); - } - else - _gBestScore = Optimum(_gBestScore, population[i]); - - for(int j = 0; j < _chromlen; ++j) - _current_position[i][j] += (float) (Configuration.NextGaussian() * 0.01 * S * (current_position[i][j] - _gBestScore[j])); - _current_position[i] = Optimum(_current_position[i], population[i]); - } - } - - private void UpdatePosition2(List population) + private void UpdateVelocities(List population) { var current_position = _current_position.ToArray(); for (int i = 0; i < _populationSize; ++i) { @@ -136,7 +77,7 @@ private void UpdatePosition2(List population) } if(changed) - _current_position[i] = Optimum(_current_position[i], population[i]); + _current_position[i] = _lf.Optimum(_current_position[i], population[i]); } } @@ -151,8 +92,8 @@ protected override void Reform() protected override List Replacement(List population) { - UpdatePosition1(population); - UpdatePosition2(population); + _gBest = _lf.UpdateVelocities(population, _populationSize, _current_position, _gBest); + UpdateVelocities(population); for (int i = 0; i < _populationSize; ++i) { var chromosome = _prototype.MakeEmptyFromPrototype(); diff --git a/GaSchedule.Algorithm/Dlba.cs b/GaSchedule.Algorithm/Dlba.cs new file mode 100644 index 0000000..049baf4 --- /dev/null +++ b/GaSchedule.Algorithm/Dlba.cs @@ -0,0 +1,222 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using GaSchedule.Model; + +/* + * Xie, Jian & Chen, Huan. (2013). + * A Novel Bat Algorithm Based on Differential Operator and Lévy Flights Trajectory. + * Computational intelligence and neuroscience. 2013. 453812. 10.1155/2013/453812. + * Copyright (c) 2024 Miller Cy Chan + */ + +namespace GaSchedule.Algorithm +{ + public class Dlba : NsgaIII where T : Chromosome + { + private int _currentGeneration, _max_iterations = 5000; + + private int _chromlen, _minValue = 0; + + private double _alpha, _pa; + + private float[] _frequency, _rate, _loudness; + + private float[] _gBest = null; + private float[][] _position = null; + private float[][] _velocity = null; + + private List _maxValues; + private LévyFlights _lf; + + // Initializes Bat algorithm + public Dlba(T prototype, int numberOfCrossoverPoints = 2, int mutationSize = 2, float crossoverProbability = 80, float mutationProbability = 3) : base(prototype, numberOfCrossoverPoints, mutationSize, crossoverProbability, mutationProbability) + { + _alpha = 0.9; + _pa = .25; + } + + static E[][] CreateArray(int rows, int cols) + { + E[][] array = new E[rows][]; + for (int i = 0; i < array.GetLength(0); i++) + array[i] = new E[cols]; + + return array; + } + + protected override void Initialize(List population) + { + _maxValues = new(); + _prototype.MakeEmptyFromPrototype(_maxValues); + + for (int i = 0; i < _populationSize; ++i) { + List positions = new(); + + // initialize new population with chromosomes randomly built using prototype + population.Add(_prototype.MakeNewFromPrototype(positions)); + + if(i < 1) { + _chromlen = positions.Count; + _frequency = new float[_chromlen]; + _rate = new float[_populationSize]; + _loudness = new float[_populationSize]; + _position = CreateArray(_populationSize, _chromlen); + _velocity = CreateArray(_populationSize, _chromlen); + _lf = new LévyFlights(_chromlen); + } + + _rate[i] = (float) Configuration.Random(); + _loudness[i] = (float) Configuration.Random() + 1; + } + } + + protected override void Reform() + { + Configuration.Seed(); + if (_crossoverProbability < 95) + _crossoverProbability += 1.0f; + else if (_pa < .5) + _pa += .01; + } + + private void UpdateVelocities(List population) + { + var mean = _loudness.Average(); + + var globalBest = _prototype.MakeEmptyFromPrototype(); + globalBest.UpdatePositions(_gBest); + var localBest = _prototype.MakeNewFromPrototype(); + + for (int i = 0; i < _populationSize; ++i) { + var beta = (float) Configuration.Random(); + var rand = Configuration.Random(); + var n = Configuration.Rand(-1.0, 1.0); + + int dim = _velocity[i].Length; + for(int j = 0; j < dim; ++j) { + _frequency[j] = ((_maxValues[j] - _minValue) * _currentGeneration / (float) n + _minValue) * beta; + _velocity[i][j] += (_position[i][j] - _gBest[j]) * _frequency[j]; + + if (rand > _rate[i]) { + _position[i][j] += _velocity[i][j]; + if(_position[i][j] > _maxValues[j]) { + _position[i][j] = _maxValues[j]; + _velocity[i][j] = _minValue; + } + else if(_position[i][j] < _minValue) + _position[i][j] = _velocity[i][j] = _minValue; + } + } + + var localTemp = _prototype.MakeEmptyFromPrototype(); + localTemp.UpdatePositions(_position[i]); + if (localTemp.Dominates(localBest)) + localBest = localTemp; + } + + for (int i = 0; i < _populationSize; ++i) { + var positionTemp = _position.ToArray(); + var rand = Configuration.Random(); + if (rand < _loudness[i]) { + var n = Configuration.Rand(-1.0, 1.0); + int dim = _velocity[i].Length; + for(int j = 0; j < dim; ++j) { + positionTemp[i][j] = _gBest[j] + (float) n * mean; + if(positionTemp[i][j] > _maxValues[j]) { + positionTemp[i][j] = _maxValues[j]; + _velocity[i][j] = _minValue; + } + else if(positionTemp[i][j] < _minValue) + positionTemp[i][j] = _velocity[i][j] = _minValue; + } + + if (globalBest.Dominates(localBest)) { + _position[i] = positionTemp[i]; + _rate[i] *= (float) Math.Pow(_currentGeneration / n, 3); + _loudness[i] *= (float) _alpha; + } + } + + _position[i] = _lf.Optimum(_position[i], population[i]); + } + } + + protected override List Replacement(List population) + { + _gBest = _lf.UpdateVelocities(population, _populationSize, _position, _gBest); + UpdateVelocities(population); + + for (int i = 0; i < _populationSize; ++i) { + var chromosome = _prototype.MakeEmptyFromPrototype(); + chromosome.UpdatePositions(_position[i]); + population[i] = chromosome; + } + + return base.Replacement(population); + } + + // Starts and executes algorithm + public override void Run(int maxRepeat = 9999, double minFitness = 0.999) + { + if (_prototype == null) + return; + + var pop = new List[2]; + pop[0] = new List(); + Initialize(pop[0]); + + // Current generation + int currentGeneration = 0; + int bestNotEnhance = 0; + double lastBestFit = 0.0; + + int cur = 0, next = 1; + while(currentGeneration < _max_iterations) + { + var best = Result; + if (currentGeneration > 0) + { + var status = string.Format("\rFitness: {0:F6}\t Generation: {1}", best.Fitness, currentGeneration); + Console.Write(status); + + // algorithm has reached criteria? + if (best.Fitness > minFitness) + break; + + var difference = Math.Abs(best.Fitness - lastBestFit); + if (difference <= 1e-6) + ++bestNotEnhance; + else { + lastBestFit = best.Fitness; + bestNotEnhance = 0; + } + + if (bestNotEnhance > (maxRepeat / 100)) + Reform(); + } + + /******************* crossover *****************/ + var offspring = Crossing(pop[cur]); + + /******************* mutation *****************/ + foreach (var child in offspring) + child.Mutation(_mutationSize, _mutationProbability); + + pop[cur].AddRange(offspring); + + /******************* replacement *****************/ + pop[next] = Replacement(pop[cur]); + _best = pop[next][0].Dominates(pop[cur][0]) ? pop[next][0] : pop[cur][0]; + + (cur, next) = (next, cur); + ++currentGeneration; + } + } + + public override string ToString() + { + return "Bat algorithm with differential operator and Levy flights trajectory (DLBA)"; + } + } +} diff --git "a/GaSchedule.Algorithm/L\303\251vyFlights.cs" "b/GaSchedule.Algorithm/L\303\251vyFlights.cs" new file mode 100644 index 0000000..0fbb4df --- /dev/null +++ "b/GaSchedule.Algorithm/L\303\251vyFlights.cs" @@ -0,0 +1,81 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using GaSchedule.Model; + +namespace GaSchedule.Algorithm +{ + internal sealed class LévyFlights where T : Chromosome { + + private int _chromlen; + private double _beta, _σu, _σv; + + internal LévyFlights(int chromlen) + { + _chromlen = chromlen; + + _beta = 1.5; + var num = Gamma(1 + _beta) * Math.Sin(Math.PI * _beta / 2); + var den = Gamma((1 + _beta) / 2) * _beta * Math.Pow(2, (_beta - 1) / 2); + _σu = Math.Pow(num / den, 1 / _beta); + _σv = 1; + } + + private static double Gamma(double z) + { + if (z < 0.5) + return Math.PI / Math.Sin(Math.PI * z) / Gamma(1.0 - z); + + // Lanczos approximation g=5, n=7 + var coef = new double[7] { 1.000000000190015, 76.18009172947146, -86.50532032941677, + 24.01409824083091, -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5 }; + + var zz = z - 1.0; + var b = zz + 5.5; // g + 0.5 + var sum = coef[0]; + for (int i = 1; i < coef.Length; ++i) + sum += coef[i] / (zz + i); + + var LogSqrtTwoPi = 0.91893853320467274178; + return Math.Exp(LogSqrtTwoPi + Math.Log(sum) - b + Math.Log(b) * (zz + 0.5)); + } + + internal float[] Optimum(float[] localVal, T chromosome) + { + var localBest = chromosome.MakeEmptyFromPrototype(); + localBest.UpdatePositions(localVal); + + if(localBest.Dominates(chromosome)) { + chromosome.UpdatePositions(localVal); + return localVal; + } + + var positions = new float[_chromlen]; + chromosome.ExtractPositions(positions); + return positions; + } + + internal float[] UpdateVelocities(List population, int populationSize, float[][] currentPosition, float[] gBest) + { + var current_position = currentPosition.ToArray(); + for(int i = 0; i < populationSize; ++i) { + var u = Configuration.NextGaussian() * _σu; + var v = Configuration.NextGaussian() * _σv; + var S = u / Math.Pow(Math.Abs(v), 1 / _beta); + + if(gBest == null) { + gBest = new float[_chromlen]; + population[i].ExtractPositions(gBest); + } + else + gBest = Optimum(gBest, population[i]); + + for(int j = 0; j < _chromlen; ++j) + currentPosition[i][j] += (float) (Configuration.NextGaussian() * 0.01 * S * (current_position[i][j] - gBest[j])); + + currentPosition[i] = Optimum(currentPosition[i], population[i]); + } + return gBest; + } + } +} diff --git a/GaSchedule.Algorithm/NsgaIII.cs b/GaSchedule.Algorithm/NsgaIII.cs index 605cdc2..6a735aa 100644 --- a/GaSchedule.Algorithm/NsgaIII.cs +++ b/GaSchedule.Algorithm/NsgaIII.cs @@ -499,14 +499,14 @@ protected virtual void Reform() { Configuration.Seed(); if (_crossoverProbability < 95) - _crossoverProbability += 1.0f; + _crossoverProbability += 1.0f; else if (_mutationProbability < 30) - _mutationProbability += 1.0f; + _mutationProbability += 1.0f; } protected virtual List Replacement(List population) { - var rps = new List(); + var rps = new List(); ReferencePoint.GenerateReferencePoints(rps, _prototype.Objectives.Length, _objDivision); return Selection(population, rps); }