Add GSA(A Gravitational Search Algorithm);

After the comparision between this code and GSA code(https://github.com/ravexina/GSA/blob/master/GSA.py), I found that the avg is 0.05 vs. 0.04 and the std is 0.04 vs 0.05 after 100 computation.

Author: LucXiong

GSA.py

@@ -0,0 +1,118 @@
+# ！usr/bin/env python
+# -*- coding: utf-8 -*-
+# Time : 2022/10/19 20:32
+# @Author : LucXiong
+# @Project : Model
+# @File : GSA.py
+
+'''
+Ref:https://github.com/ravexina/GSA/blob/master/GSA.py
+Ref:Rashedi E., Nezamabadi-Pour H., Saryazdi S. GSA: A Gravitational Search Algorithm[J]. Information Sciences, 2009, 179(13): 2232-48.
+'''
+
+import numpy as np
+import matplotlib.pyplot as plt
+import math
+import test_function
+
+class GSA():
+    def __init__(self, func, n_dim=None, pop=40, max_iter=150, lb=-1e5, ub=1e5, alpha=0.1, G=0.9):
+        self.func = func
+        self.alpha = alpha
+        self.G = G
+        self.pop = pop  # number of particles
+        self.n_dim = n_dim  # dimension of particles, which is the number of variables of func
+        self.max_iter = max_iter  # max iter
+
+        self.lb, self.ub = np.array(lb) * np.ones(self.n_dim), np.array(ub) * np.ones(self.n_dim)
+        assert self.n_dim == len(self.lb) == len(self.ub), 'dim == len(lb) == len(ub) is not True'
+        assert np.all(self.ub > self.lb), 'upper-bound must be greater than lower-bound'
+
+        self.X = np.random.uniform(low=self.lb, high=self.ub, size=(self.pop, self.n_dim))
+        v_high = (self.ub - self.lb) # 速度设置为区间长度的一半
+        self.V = np.random.uniform(low=-v_high, high=v_high, size=(self.pop, self.n_dim))  # speed of particles
+        self.Y = [self.func(self.X[i]) for i in range(self.pop)]  # y = f(x) for all particles
+        self.q = [1 for i in range(self.pop)]
+        self.M = [1 for i in range(self.pop)]
+        self.f = [[0 for j in range(self.n_dim)] for i in range(self.pop)]
+        self.a = [[0 for j in range(self.n_dim)] for i in range(self.pop)]
+
+
+        self.pbest_x = self.X.copy()  # personal best location of every particle in history
+        self.pbest_y = [np.inf for i in range(self.pop)]  # best image of every particle in history
+        self.gbest_x = self.pbest_x.mean(axis=0).reshape(1, -1)  # global best location for all particles
+        self.gbest_y = np.inf  # global best y for all particles
+        self.gbest_y_hist = []  # gbest_y of every iteration
+        self.update_gbest()
+
+    def cal_q_M(self):
+        best = np.min(self.Y)
+        worst = np.max(self.Y)
+        self.q = (self.Y - worst) / best - worst
+        self.M = self.q / sum(self.q)
+
+    def cal_f(self):
+        for i in range(self.pop):
+            f = None
+            for j in range(self.pop):
+                if j != i:
+                    dividend = float(self.M[i] * self.M[j])
+                    temp = self.X[i] - self.X[j]
+                    sum_temp = [k**2 for k in temp]
+                    divisor = math.sqrt(sum(sum_temp)) + np.finfo('float').eps
+                    if f is None:
+                        f = self.G * (dividend / divisor) * (self.X[j] - self.X[i])
+                    else:
+                        f = f + self.G * (dividend / divisor) * (self.X[j] - self.X[i])
+
+            self.f[i] = np.random.uniform(0, 1) * f
+
+    def update_gbest(self):
+        idx_min = self.Y.index(min(self.Y))
+        if self.gbest_y > self.Y[idx_min]:
+            self.gbest_x = self.X[idx_min, :].copy()
+            self.gbest_y = self.Y[idx_min]
+
+    def run(self):
+        for iteration in range(self.max_iter):
+
+            self.Y = [self.func(self.X[i]) for i in range(self.pop)]
+            self.cal_q_M()
+            self.G = self.G * np.e ** (- self.alpha * (iteration / self.max_iter))
+            self.cal_f()
+            self.a = [self.f[i]/self.M[i] for i in range(self.pop)]
+            self.V = (np.random.uniform(0, 1) * self.V) + self.a
+            self.update_gbest()
+            self.X = self.X + self.V
+            self.gbest_y_hist.append(self.gbest_y)
+            # print(iteration, self.gbest_x, self.gbest_y)
+        self.best_x, self.best_y = self.gbest_x, self.gbest_y
+        return self.best_x, self.best_y
+
+
+def demo_func(args):
+
+    x, y = args[0], args[1]
+    a = 1
+    b = 100
+    return (a - x) ** 2 + b * (y - x ** 2) ** 2
+
+
+if __name__ == '__main__':
+    n_dim = 2
+    lb = [0 for i in range(n_dim)]
+    ub = [1 for i in range(n_dim)]
+    # demo_func = test_function.fu2
+    pop_size = 20
+    max_iter = 100
+    res = []
+    for i in range(100):
+        pso = GSA(func=demo_func, n_dim=n_dim, pop=pop_size, max_iter=max_iter, lb=lb, ub=ub)
+        best_x, bext_y = pso.run()
+        print(f'{i}: {demo_func(pso.gbest_x)}\t{pso.gbest_x}')
+        res.append(bext_y)
+    print(sum(res)/len(res))
+    print(np.std(res))
+    # plt.plot(pso.gbest_y_hist)
+    #
+    # plt.show()