Potential Incorrect Fitness issue ?

[vsedov]

Im rather new to optimisation so forgive me if i sound confused. But ive been trying out gwo with ackley function

population = 50 dimension = 30 position_min = -15 position_max = 30 generation = 500 fitness_error = 1e-18

Issue is, running on all other grey wolf algos ive tested i average around

Starting GWO algorithm Number of wolves (N): 50 Dimension (dim): 30 Max cycles (iterations): 500 Search space bounds: -15 to 30

7.549516567451064e-15 For best Fitness.

-15 or -16, while this variation, i can barely get -02.

[vsedov]

Ive gone through the code base, and the gwo algo seems alright to me, i couldnt find any issues there; hence my confusion as to why this might be occurring

[matteocacciola]

Hello @vsedov , can you please provide the code you are trying to implement, with some more details about the possible issue you are facing?

[vsedov]

Hi,

Thanks for getting back to me;

Problem case;

For ackley, im currently working on GWO; and working on modifying it; Just to see the average baseline of gwo i wanted to see what it would look like for

population = 50
dimension = 30
position_min = -15
position_max = 30
generation = 500
fitness_error = 1e-18```Following settings; i noticed you updated your code base, but the key factor here is dimensions being 30. 

Provided test code i found / made

I did this to compare as i had some suspicion that my current results were not of what i was expecting.

import math
import random

import numpy as np

def vec(dim):
    return [0.0 for i in range(dim)]

class FitnessModel:
    """
    fitness_function    : Target Model (unimodal, multimodal, fix_dim_multimodal or composite)
    _range              : Range of fitness function
    """

    def __init__(self, fitness_function, _range):
        self.eval = fitness_function
        self.xmin = _range[0]
        self.xmax = _range[1]

class Wolf:
    """
    fitness : FitnessModel object function to evaluate fitness
    dim     : dimension of the wolf position vector
    """

    def __init__(self, fitness, dim, seed):
        # self.rnd = random.Random(seed)
        self.rnd = np.random.RandomState(seed)

        self.position = vec(dim)

        for i in range(dim):
            self.position[i] = (
                fitness.xmax - fitness.xmin
            ) * self.rnd.random() + fitness.xmin

        self.fitness = fitness.eval(self.position)

class GWO:
    """
    fitness : fitness model
    max_iter: maximum number of iteration
    N       : number of wolves/search agents.
    dim     : dimension of the position vector
    """

    def __init__(self, fitness, max_iter, N, dim, seed):
        self.fitness = fitness
        self.max_iter = max_iter
        self.num_wolves = N
        self.dim = dim
        self.rnd = np.random.RandomState(seed)
        self.X = [Wolf(fitness, dim, i) for i in range(num_wolves)]
        self.alpha, self.beta, self.delta, *self.omega = self.sortAgents()
        self.a = 2

    def sortAgents(self):
        self.X = sorted(self.X, key=lambda wolf: wolf.fitness)
        return self.X

    def init_vec_A(self):
        # A = vec(3)
        # for i in range(3):
        #     # A[i] = self.a * (2 * self.rnd.random() - 1)
        #     A[i] = self.a * (2 * np.random.random(3) -1)
        # return A
        return self.a * (2 * np.random.random(3) - 1)

    def init_vec_C(self):
        C = vec(3)
        for i in range(3):
            C[i] = 2 * self.rnd.random()
        return C

    def optimize(self):

        t = 0
        while t < self.max_iter:
            # print(
            #     "Iter = "
            #     + str(t)
            #     + " alpha position = "
            #     + str(self.alpha.position)
            #     # + " alpha fitness = "  self.alpha.fitness
            # )

            a = 2 * (1 - t / self.max_iter)

            for i in range(self.num_wolves):
                A1, A2, A3 = self.init_vec_A()
                C1, C2, C3 = self.init_vec_C()

                D_alpha = vec(self.dim)
                D_beta = vec(self.dim)
                D_delta = vec(self.dim)

                for j in range(self.dim):
                    D_alpha[j] = np.abs(
                        C1 * self.alpha.position[j] - self.X[i].position[j]
                    )
                    D_beta[j] = np.abs(
                        C2 * self.beta.position[j] - self.X[i].position[j]
                    )
                    D_delta[j] = np.abs(
                        C3 * self.delta.position[j] - self.X[i].position[j]
                    )

                X1 = vec(self.dim)
                X2 = vec(self.dim)
                X3 = vec(self.dim)

                for j in range(self.dim):
                    X1[j] = self.alpha.position[j] - A1 * D_alpha[j]
                    X2[j] = self.beta.position[j] - A2 * D_beta[j]
                    X3[j] = self.delta.position[j] - A3 * D_delta[j]

                Xnew = vec(self.dim)

                for j in range(self.dim):
                    Xnew[j] = (X1[j] + X2[j] + X3[j]) / 3

                fnew = self.fitness.eval(Xnew)
                if fnew < self.X[i].fitness:
                    self.X[i].position = Xnew
                    self.X[i].fitness = fnew

            self.alpha, self.beta, self.delta, *self.omega = self.sortAgents()

            t += 1
        return self.alpha.position
def aclkey(x):
    A = 20
    B = 0.2
    C = 2 * np.pi
    d = len(x)  # Correct handling of dimensionality
    return (
        -A * np.exp(-B * np.sqrt(sum([xi**2 for xi in x]) / d))
        - np.exp(sum([np.cos(C * xi) for xi in x]) / d)
        + A
        + np.exp(1)
    )

""" ----------------------------------------------------------- """

runs = []

dim = 30
fitness = FitnessModel(fitness_function=aclkey, _range=[-15, 30])

num_wolves = 50
max_iter = 500

print("num_wolves = " + str(num_wolves))
print("max_iter = " + str(max_iter))
print("\nStarting GWO algorithm\n")

optimiser = GWO(fitness, max_iter, num_wolves, dim, 0)
optimum_position = optimiser.optimize()

print("\nGWO completed\n")
print("\nBest solution found:")
print([optimum_position[k] for k in range(dim)])

err = fitness.eval(optimum_position)
print("\n\n\n\n\n")
print("fitness of best solution => ", err)
runs.append(err)

print("\n\n\n\n\n")

I was expecting somewhere around 5.28E-12 to -15 right as a baseline.

With Dimensions of 30

There is a potential chance im doing something wrong; i also tested this idea on the following code too

import copy  # array-copying convenience
import math  # cos() for Rastrigin
import random
import sys  # max float

# -------fitness functions---------

# rastrigin function
#
def fitness_rastrigin(position):
    fitness_value = 0.0
    for i in range(len(position)):
        xi = position[i]
        fitness_value += (xi * xi) - (10 * math.cos(2 * math.pi * xi)) + 10
    return fitness_value

# sphere function
def fitness_sphere(position):
    fitness_value = 0.0
    for i in range(len(position)):
        xi = position[i]
        fitness_value += xi * xi
    return fitness_value

# -------------------------

# wolf class
class wolf:
    def __init__(self, fitness, dim, minx, maxx, seed):
        self.rnd = random.Random(seed)
        self.position = [0.0 for i in range(dim)]
        for i in range(dim):
            self.position[i] = (maxx - minx) * self.rnd.random() + minx
        self.fitness = fitness(self.position)  # curr fitness

def gwo(fitness, max_iter, n, dim, minx, maxx):
    rnd = random.Random(0)

    # create n random wolves
    population = [wolf(fitness, dim, minx, maxx, i) for i in range(n)]

    # On the basis of fitness values of wolves
    # sort the population in asc order
    population = sorted(population, key=lambda temp: temp.fitness)

    # best 3 solutions will be called as
    # alpha, beta and gaama
    alpha_wolf, beta_wolf, gamma_wolf = copy.copy(population[:3])

    # main loop of gwo
    Iter = 0
    while Iter < max_iter:

        # after every 10 iterations
        # print iteration number and best fitness value so far
        if Iter % 10 == 0 and Iter > 1:
            print(
                "Iter = "
                + str(Iter)
                + " best fitness = %.3f" % alpha_wolf.fitness
            )

        # linearly decreased from 2 to 0
        a = 2 * (1 - Iter / max_iter)

        # updating each population member with the help of best three members
        for i in range(n):
            A1, A2, A3 = (
                a * (2 * rnd.random() - 1),
                a * (2 * rnd.random() - 1),
                a * (2 * rnd.random() - 1),
            )
            C1, C2, C3 = (
                2 * rnd.random(),
                2 * rnd.random(),
                2 * rnd.random(),
            )

            X1 = [0.0 for i in range(dim)]
            X2 = [0.0 for i in range(dim)]
            X3 = [0.0 for i in range(dim)]
            Xnew = [0.0 for i in range(dim)]
            for j in range(dim):
                X1[j] = alpha_wolf.position[j] - A1 * abs(
                    C1 * alpha_wolf.position[j] - population[i].position[j]
                )
                X2[j] = beta_wolf.position[j] - A2 * abs(
                    C2 * beta_wolf.position[j] - population[i].position[j]
                )
                X3[j] = gamma_wolf.position[j] - A3 * abs(
                    C3 * gamma_wolf.position[j] - population[i].position[j]
                )
                Xnew[j] += X1[j] + X2[j] + X3[j]

            for j in range(dim):
                Xnew[j] /= 3.0

            # fitness calculation of new solution
            fnew = fitness(Xnew)

            # greedy selection
            if fnew < population[i].fitness:
                population[i].position = Xnew
                population[i].fitness = fnew

        # On the basis of fitness values of wolves
        # sort the population in asc order
        population = sorted(population, key=lambda temp: temp.fitness)

        # best 3 solutions will be called as
        # alpha, beta and gaama
        alpha_wolf, beta_wolf, gamma_wolf = copy.copy(population[:3])

        Iter += 1
    # end-while

    # returning the best solution
    return alpha_wolf.position

# ----------------------------

# Driver code for rastrigin function

print("\nBegin grey wolf optimization on rastrigin function\n")
dim = 30

def ackley(x):
    a = 20
    b = 0.2
    c = 2 * math.pi
    s1 = sum([i**2 for i in x])
    s2 = sum([math.cos(c * i) for i in x])
    return (
        -a * math.exp(-b * math.sqrt(s1 / len(x)))
        - math.exp(s2 / len(x))
        + a
        + math.exp(1)
    )

fitness = ackley
num_particles = 50
max_iter = 500

print("Setting num_particles = " + str(num_particles))
print("Setting max_iter = " + str(max_iter))
print("\nStarting GWO algorithm\n")

best_position = gwo(fitness, max_iter, num_particles, dim, -15.0, 30.0)

print("\nGWO completed\n")
print("\nBest solution found:")
print([best_position[k] for k in range(dim)])
err = fitness(best_position)
print("fitness of best solution =", err)

print("\nEnd GWO for rastrigin\n")

import math
import random

import numpy as np

class GreyWolf:
    def __init__(self, position, fitness):
        self.position = np.array(position)
        self.fitness = fitness

class GreyWolfOptimization:
    def __init__(self, fitness_function, config=None, debug=False):
        self.fitness_function = fitness_function
        self.config = config
        self.debug = debug
        self.population = []
        self.alpha_wolf = None
        self.beta_wolf = None
        self.delta_wolf = None
        self.initialize_wolves()

    def initialize_wolves(self):
        for _ in range(self.config.N):
            position = np.random.uniform(
                self.config.lower_bound,
                self.config.upper_bound,
                self.config.dim,
            )
            fitness = self.fitness_function(position)
            wolf = GreyWolf(position, fitness)
            self.population.append(wolf)
        self.update_leading_wolves()

    def update_leading_wolves(self):
        self.population.sort(key=lambda wolf: wolf.fitness)
        self.alpha_wolf, self.beta_wolf, self.delta_wolf = self.population[:3]

    def evolve(self):
        a = 2 - self.current_cycle * (
            2 / self.config.max_cycles
        )  # Decrease linearly from 2 to 0
        for wolf in self.population:
            for i in range(self.config.dim):
                A1, A2, A3 = (
                    2 * a * np.random.rand(3) - a
                )  # Coefficients for alpha, beta, and delta
                C1, C2, C3 = 2 * np.random.rand(3)

                D_alpha = np.abs(
                    C1 * self.alpha_wolf.position[i] - wolf.position[i]
                )
                D_beta = np.abs(
                    C2 * self.beta_wolf.position[i] - wolf.position[i]
                )
                D_delta = np.abs(
                    C3 * self.delta_wolf.position[i] - wolf.position[i]
                )

                X1 = self.alpha_wolf.position[i] - A1 * D_alpha
                X2 = self.beta_wolf.position[i] - A2 * D_beta
                X3 = self.delta_wolf.position[i] - A3 * D_delta

                wolf.position[i] = (X1 + X2 + X3) / 3

            wolf.fitness = self.fitness_function(wolf.position)

    def optimize(self):
        self.current_cycle = 0
        if self.debug:
            print("Starting GWO algorithm")
            print(f"Number of wolves (N): {self.config.N}")
            print(f"Dimension (dim): {self.config.dim}")
            print(f"Max cycles (iterations): {self.config.max_cycles}")
            print(
                f"Search space bounds: {self.config.lower_bound} to {self.config.upper_bound}\n"
            )
            print("Starting optimization with GWO...")
        while self.current_cycle < self.config.max_cycles:
            self.evolve()
            self.update_leading_wolves()
            self.current_cycle += 1
            if self.debug:
                print(
                    f"Cycle: {self.current_cycle}, Alpha Fitness: {self.alpha_wolf.fitness}"
                )

        if self.debug:
            print("\nGWO completed")
            print("Best solution found:", self.alpha_wolf.position)
            print("Fitness of the best solution:", self.alpha_wolf.fitness)
        return self.alpha_wolf.position, self.alpha_wolf.fitness

# Definition of the Ackley function
def ackley_function(x):
    A = 20
    B = 0.2
    C = 2 * np.pi
    d = len(x)  # Correct handling of dimensionality
    return (
        -A * np.exp(-B * np.sqrt(sum([xi**2 for xi in x]) / d))
        - np.exp(sum([np.cos(C * xi) for xi in x]) / d)
        + A
        + np.exp(1)
    )

class GreyWolfOptimizationConfig:
    def __init__(
        self,
        N=50,
        dim=30,
        max_cycles=500,
        lower_bound=-32.768,
        upper_bound=32.768,
    ):
        self.N = N  # Number of wolves
        self.dim = dim  # Dimension of the problem
        self.max_cycles = max_cycles  # Number of iterations
        self.lower_bound = lower_bound  # Lower bound of the search space
        self.upper_bound = upper_bound  # Upper bound of the search space

config = GreyWolfOptimizationConfig(
    # dim=30, max_cycles=500, lower_bound=-15, upper_bound=30,
    N=50,
    dim=30,
    max_cycles=500,
    lower_bound=-15,
    upper_bound=30,
)

gwo = GreyWolfOptimization(
    fitness_function=ackley_function, config=config, debug=True
)

best_position, best_fitness = gwo.optimize()

# Print  best solution found
# print(f"Best Position: {best_position}")
for i in range(len(best_position)):
    print(f"x{i} = {best_position[i]}")

print(f"Best Fitness: {best_fitness}")

As for what i was doing with the current code base with pyevo :

from demos.functions.ackley import (
    fitness_error,
    generation,
    name,
    population,
    position_max,
    position_min,
    task,
)
from pyvolutionary import GreyWolfOptimization, GreyWolfOptimizationConfig
from pyvolutionary.utils import animate

configuration = GreyWolfOptimizationConfig(
    population_size=population,
    fitness_error=fitness_error,
    max_cycles=generation,
)

optimization_result = GreyWolfOptimization(configuration, True).optimize(task)

print(data.best_solution.fitness)
print(data.best_solution.cost)

#########
import numpy as np
from pyvolutionary import Task, ContinuousVariable

class Ackley(Task):
    def objective_function(self, x: list[float]) -> float:
        A = 20
        B = 0.2
        C = 2 * np.pi
        return -A * np.exp(-B * np.sqrt(sum([xi ** 2 for xi in x]) / dimension)) - np.exp(
            sum([np.cos(C * xi) for xi in x]) / dimension
        ) + A + np.exp(1)

population =30
dimension = 30
position_min = -15.0
position_max = 30.0
generation = 500

task = Ackley(
    variables=[ContinuousVariable(
        name=f"x{i}", lower_bound=position_min, upper_bound=position_max
    ) for i in range(dimension)],
)
name = "ackley"

Fitness error was removed as i had modified the code to where i wanted to run the full iterations.

I also noticed that you had changed the code base, so will need to see if i can replicate similar results somehow.

Extra information

I suspect its just something to do with how ever the hell im defining the dims;

Just in case, how would we be able to set dimensions right now ? if say i wanted to have a test function with 30 dim on ackley for example.

that being said, thanks once again for getting back to me and also great work on this project, really impressive. Il mess around a bit more as im almost certain im being stupid here but thought id ask just in case. Il also try messing around with opfunu perhaps as well. If you have any idea what im doing wrong please do let me know, would be great to know.

Guess im just confused why i cant seem to replicate results or even near to it.

Thanks once again

[matteocacciola]

Hello, I see a lot of code you implemented by yourself. I will be happy to help you with the code using pyVolutionary and not working. Can you please provide that code only?

[vsedov]

Hi,

Yeah so the code i gave is just some test / benchmarks i was trying to replicate with pyVolutionary if that makes sense, running those code bases; all provided similar results.

What i have right now with pyVolutionary is the following:

Ackley function

import numpy as np
from pyvolutionary import ContinuousMultiVariable, Task

class Ackley(Task):
    def objective_function(self, x: list[float]) -> float:
        A = 20
        B = 0.2
        C = 2 * np.pi
        dimension = len(x)
        return (
            -A * np.exp(-B * np.sqrt(sum([xi**2 for xi in x]) / dimension))
            - np.exp(sum([np.cos(C * xi) for xi in x]) / dimension)
            + A
            + np.exp(1)
        )

population = 100
generation = 500
fitness_error = 1e-100
task = Ackley(
    variables=[
        ContinuousMultiVariable(
            name="x",
            lower_bounds=[-15] * 30,
            upper_bounds=[30] * 30,
        )
    ],
)
name = "ackley"

How im running it

from demos.functions.ackley import generation, name, population, task
from pyvolutionary import GreyWolfOptimization, GreyWolfOptimizationConfig
from pyvolutionary.utils import animate

configuration = GreyWolfOptimizationConfig(
    population_size=population,
    fitness_error=1e-100,
    max_cycles=generation,
)

data = GreyWolfOptimization(configuration, True).optimize(task)

print(data.best_solution.fitness)
print(data.best_solution.cost)

Pretty much stemming from the demo .

Potential Issue Within Grey wolf algo

My current implementation; With all my wolves is the following: This is something i just modified, and the implmentation used in the previous examples too, and somethign i double checked with the paper as well.

    def optimization_step(self):
        def evolve(wolf: GreyWolf) -> GreyWolf:
            pos = np.array(wolf.position)
            # linearly decreased from 2 to 0
    ################ FIX
            a = 2 * (1 - self._current_cycle / self._config.max_cycles)
            a1, a2, a3 = (a * (2 * np.random.random(3) - 1)).tolist()
            c1, c2, c3 = (2 * np.random.random(3)).tolist()
    ################ FIX
            x1 = np.array(self.__alpha_wolf.position) - a1 * np.abs(
                c1 * np.array(self.__alpha_wolf.position) - pos
            )
            x2 = np.array(self.__beta_wolf.position) - a2 * np.abs(
                c2 * np.array(self.__beta_wolf.position) - pos
            )
            x3 = np.array(self.__gamma_wolf.position) - a3 * np.abs(
                c3 * np.array(self.__gamma_wolf.position) - pos
            )
            # greedy selection
            return self._greedy_select_agent(
                wolf,
                GreyWolf(**self._init_agent((x1 + x2 + x3) / 3).model_dump()),
            )

        # updating each population member with the help of best three members
        self._population = [evolve(wolf) for wolf in self._population]

        # best 3 solutions will be called as alpha, beta and gaama
        self.__alpha_wolf, self.__beta_wolf, self.__gamma_wolf = best_agents(
            self._population, 3
        )

Could you please verify with your current implementation of grey wolf and the current modification i made is correct?

just to be sure. I beleive there is an error in teh grey wolf optimisation_step function.

Original Function :

    def optimization_step(self):
        def evolve(wolf: GreyWolf) -> GreyWolf:
            pos = np.array(wolf.position)
            x1 = np.array(self.__alpha_wolf.position) - a1 * np.abs(
                c1 * np.array(self.__alpha_wolf.position) - pos
            )
            x2 = np.array(self.__beta_wolf.position) - a2 * np.abs(
                c2 * np.array(self.__beta_wolf.position) - pos
            )
            x3 = np.array(self.__gamma_wolf.position) - a3 * np.abs(
                c3 * np.array(self.__gamma_wolf.position) - pos
            )
            # greedy selection
            return self._greedy_select_agent(
                wolf,
                GreyWolf(**self._init_agent((x1 + x2 + x3) / 3).model_dump()),
            )
    ################ Original
        # linearly decreased from 2 to 0
        a = 2 * (1 - self._current_cycle / self._config.max_cycles)
        a1, a2, a3 = (a * (2 * np.random.random(3) - 1)).tolist()
        c1, c2, c3 = (2 * np.random.random(3)).tolist()
    ################ Original
        # updating each population member with the help of best three members
        self._population = [evolve(wolf) for wolf in self._population]

        # best 3 solutions will be called as alpha, beta and gaama
        self.__alpha_wolf, self.__beta_wolf, self.__gamma_wolf = best_agents(
            self._population, 3
        )

Hope this helps; but after changing that, the issue resolved it self and i was getting the results i needed.

It seems your current implementation the coeficients are calculated once per optimisation the for the entire population, while the modified one i just showed is for each wolf, which i think is the correct implementation. Would you mind verifying this for me ?

And just in case, you know when defining dimensions for functions with pyVolutionary What would be a good way of going about it ?

[vsedov]

So with the modified version, im getting the expected results;

Changed code

Maximum number of cycles reached best_solution.fitness 0.999999999999996 best_solution.cost 3.9968028886505635e-15

Previous Code

0.5286819064805379 0.8914965459231436

compared to the original grey wolf algo without any modifications.

[matteocacciola]

@vsedov this is the code I tested (from yours):

from demos.functions.ackley import generation, population, task
from pyvolutionary import GreyWolfOptimization, GreyWolfOptimizationConfig

configuration = GreyWolfOptimizationConfig(
    population_size=population,
    fitness_error=1e-100,
    max_cycles=generation,
)

data = GreyWolfOptimization(configuration, True).optimize(task)

print(data.best_solution.fitness)
print(data.best_solution.cost)

Here are the results I obtain

0.999999999999996
3.9968028886505635e-15

Obviously, they can change according to the initialization. If you run the optimization different times, you can see different results: it is due to the initial randomization of the agents. This is why these algorithms are called "meta-heuristic". They are able to provide sub-optimal results (i.e., results close to the exact one with an error below an acceptable threshold) in a reasonable time and a reasonable computational load. By the way, I see your proposal and I see that the randomization of a1, a2, a3, c1, c2, c3 parameters could be stressed for each iteration and for each agent, rather than for each iteration only. But the latter is the way reported by the original paper: https://seyedalimirjalili.com/gwo.

If you want to contribute to pyVolutionary, could you please propose a PR, so that I can check and evaluate your improvements? The analysis of a PR is easier and faster.

N.B. Side note: in case of usage of meta-heuristic algorithms, I always suggest to run different trials and make a sort of mathematical mean of the optimizer's performances. That's why I also provide the Multitask helper in my pyVolutionary.

[vsedov]

@matteocacciola Thanks for letting me know; just to be sure il go ahead and check the original implementation, and will create a pr if needed as after double checking my self i think you are correct there.

And also thanks for letting me know about the multitask, will be using that more from now on. Guess i was really stressed out from the results i was obtaining for a while, but yeah, il go ahead and see if there was a mistake put on my end.