遗传算法解决低维和高维函数最优化问题

遗传算法解决低维和高维函数最优化问题

import numpy as np
from random import random


class EABi:
    def __init__(self, n_var, lb, ub, ff, wf, populations=50, iterations=2000, eps=1e-4,
                 cross_method='1', p_cross=0.1, mutation_method='1', p_mutation=0.05,
                 choose_method='1', cnt_saves=4):
        self.n = n_var
        self.fitness_function = ff
        # wrapper for MIN function
        self.wrapper_function = wf
        self.populations = populations
        self.iterations = iterations
        self.eps = eps
        self.lb = lb
        self.ub = ub
        # length of each variables(in BIN)
        self.l = np.asarray(np.ceil(np.log2((self.ub - self.lb) / eps + 1)), dtype=np.int32)
        self.end_of_var = np.concatenate((np.array([0], dtype=np.int32), np.cumsum(self.l)), axis=None)
        # length of each individual
        self.len_individual = np.sum(self.l)

        self.cross_method = cross_method
        self.p_cross = p_cross
        self.mutation_method = mutation_method
        self.p_mutation = p_mutation
        self.choose_method = choose_method

        # generate pow(2, n) for speed, shape=(len, 1)
        self.pow2 = 2 ** np.arange(np.max(self.l) + 2).reshape(-1, 1)

        self.chrom, self.fitness = self.initRace()
        self.new_chrom = None

        # frequency to save chrom
        self.save_freq = int(self.iterations / cnt_saves)

    def decode(self, individual=None):
        if individual is None:
            res = np.zeros((self.populations, self.n))
            for j in range(self.n):
                res[:, j] = self.lb[j] + (self.ub[j] - self.lb[j]) / (self.pow2[self.l[j], 0] - 1) * 
                            np.dot(self.chrom[:, self.end_of_var[j]:self.end_of_var[j + 1]],
                                   self.pow2[:self.l[j]]).reshape((-1,))
        else:
            individual = individual.reshape((-1, self.len_individual))
            res = np.zeros((individual.shape[0], self.n))
            for j in range(self.n):
                res[:, j] = self.lb[j] + (self.ub[j] - self.lb[j]) / (self.pow2[self.l[j], 0] - 1) * 
                            np.dot(individual[:, self.end_of_var[j]:self.end_of_var[j + 1]],
                                   self.pow2[:self.l[j]]).reshape((-1,))
        return res

    def initRace(self):
        chrom = np.random.rand(self.populations, self.len_individual) > 0.5
        fitness = self.getFitness(chrom)
        return chrom, fitness

    def cross(self):
        if self.cross_method == '1':
            for i in range(self.populations):
                for j in range(i + 1, self.populations):
                    if random() < self.p_cross:
                        t = np.random.randint(0, self.len_individual)
                        tmp = self.new_chrom[i, t:]
                        self.new_chrom[i, t:] = self.new_chrom[j, t:]
                        self.new_chrom[j, t:] = tmp
        else:
            pass

    def mutation(self):
        if self.cross_method == '1':
            rand = np.random.rand(*self.chrom.shape)
            is_change = rand < self.p_mutation
            self.new_chrom = self.new_chrom ^ is_change
        else:
            pass

    def choose(self, keep_cnt=2):
        self.new_chrom = np.vstack((self.new_chrom, self.chrom), )
        fitness = self.getFitness()
        tmp_id = np.argsort(-fitness)
        self.new_chrom = self.new_chrom[tmp_id, :]
        fitness = fitness[tmp_id]
        # Keep the top keep_cnt-th individual
        self.chrom[0:keep_cnt, :] = self.new_chrom[0:keep_cnt, :]
        self.fitness[0:keep_cnt] = fitness[0:keep_cnt]
        # Decide the others
        if self.choose_method == '1':
            # Lun pan du
            # CAUTION: `fitness` MUST BE POSITIVE
            # adjust `self.wrapper_function` to make it
            fitness = fitness[keep_cnt:]
            self.new_chrom = self.new_chrom[keep_cnt:, ]
            p = fitness / np.sum(fitness)
            sum_p = np.cumsum(p)
            cnt = keep_cnt
            while cnt < self.populations:
                r = random()
                j = 0
                for i in range(p.shape[0]):
                    if r <= sum_p[i]:
                        j = i
                        break
                self.chrom[cnt, :] = self.new_chrom[j, :]
                self.fitness[cnt] = fitness[j]
                cnt += 1
        else:
            pass

    def getFitness(self, x=None):
        if x is None:
            de = self.decode(self.new_chrom)
            fitness = np.array(
                [self.wrapper_function(self.fitness_function(de[i, :])) for i in range(self.new_chrom.shape[0])])
        else:
            x = x.reshape((-1, self.len_individual))
            de = self.decode(x)
            fitness = np.array([self.wrapper_function(self.fitness_function(de[i, :])) for i in range(x.shape[0])])
        return fitness

    def main(self):
        best_fit = [np.max(self.fitness)]
        mean_fit = [np.mean(self.fitness)]
        saved_idx = [0]
        saved_chrom = [self.decode()]
        for i in range(1, self.iterations):
            self.new_chrom = self.chrom.copy()
            self.cross()
            self.mutation()
            self.choose()
            best_fit.append(np.max(self.fitness))
            mean_fit.append(np.mean(self.fitness))
            # print(np.mean(self.fitness), np.max(self.fitness))
            if i % self.save_freq == 0:
                saved_idx.append(i)
                saved_chrom.append(self.decode())
        saved_idx.append(self.iterations)
        saved_chrom.append(self.decode())
        return np.array(mean_fit), np.array(best_fit), saved_idx, np.array(saved_chrom)


class EAReal:
    def __init__(self, n_var, lb, ub, ff, wf, populations=50, iterations=2000, eps=1,
                 cross_method='1', p_cross=0.2, mutation_method='1', p_mutation=0.8,
                 choose_method='1', cnt_saves=4):
        self.n = n_var
        self.fitness_function = ff
        # wrapper for MIN function
        self.wrapper_function = wf
        self.populations = populations
        self.iterations = iterations
        self.eps = eps
        self.len_individual = self.n  # length of each individual
        self.lb = lb.reshape((1, -1))
        self.ub = ub.reshape((1, -1))

        self.cross_method = cross_method
        self.p_cross = p_cross
        self.mutation_method = mutation_method
        self.p_mutation = p_mutation

        self.choose_method = choose_method

        self.chrom, self.fitness = self.initRace()
        self.new_chrom = None

        # frequency to save chrom
        self.save_freq = int(self.iterations / cnt_saves)

    def decode(self):
        return self.chrom.copy()

    def initRace(self):
        chrom = np.random.rand(self.populations, self.n)
        chrom = chrom * (self.ub - self.lb) + self.lb
        return chrom, self.getFitness(chrom)

    def cross(self, lamb=0.7):
        if self.cross_method == '1':
            for i in range(self.populations):
                for j in range(i + 1, self.populations):
                    if random() < self.p_cross:
                        self.new_chrom[i, :] = lamb * self.chrom[i, :] + (1 - lamb) * self.chrom[j, :]
                        self.new_chrom[j, :] = lamb * self.chrom[j, :] + (1 - lamb) * self.chrom[i, :]
        else:
            pass

    def mutation(self):
        if self.cross_method == '1':
            rand = np.random.rand(*self.chrom.shape)
            is_change = rand < self.p_mutation
            rand = np.random.rand(*self.chrom.shape)
            is_neg = np.ones_like(is_change, dtype=int)
            is_neg[rand < 0.5] = -1
            self.new_chrom = self.new_chrom + is_change * is_neg * self.eps
        else:
            pass

    def choose(self, keep_cnt=2):
        assert np.all(self.new_chrom <= self.ub) and np.all(self.new_chrom >= self.lb)
        self.new_chrom = np.vstack((self.new_chrom, self.chrom), )
        fitness = self.getFitness()
        tmp_id = np.argsort(-fitness)
        self.new_chrom = self.new_chrom[tmp_id, :]
        fitness = fitness[tmp_id]
        # Keep the top keep_cnt-th individual
        self.chrom[0:keep_cnt, :] = self.new_chrom[0:keep_cnt, :]
        self.fitness[0:keep_cnt] = fitness[0:keep_cnt]
        # Decide the others
        if self.choose_method == '1':
            # Lun pan du
            # CAUTION: `fitness` MUST BE POSITIVE
            # adjust `self.wrapper_function` to make it
            fitness = fitness[keep_cnt:]
            self.new_chrom = self.new_chrom[keep_cnt:, ]
            p = fitness / np.sum(fitness)
            sum_p = np.cumsum(p)
            cnt = keep_cnt
            while cnt < self.populations:
                r = random()
                j = 0
                for i in range(p.shape[0]):
                    if r <= sum_p[i]:
                        j = i
                        break
                self.chrom[cnt, :] = self.new_chrom[j, :]
                self.fitness[cnt] = fitness[j]
                cnt += 1
        else:
            pass

    def getFitness(self, x=None):
        if x is None:
            fitness = np.array([self.wrapper_function(self.fitness_function(self.new_chrom[i, :])) for i in
                                range(self.new_chrom.shape[0])])
        else:
            x = x.reshape((-1, self.n))
            fitness = np.array([self.wrapper_function(self.fitness_function(x[i, :])) for i in range(x.shape[0])])
        return fitness

    def main(self):
        best_fit = [np.max(self.fitness)]
        mean_fit = [np.mean(self.fitness)]
        saved_idx = [0]
        saved_chrom = [self.decode()]
        for i in range(1, self.iterations):
            if i % 50 == 0:
                self.ub = np.max(self.chrom, axis=0, keepdims=True)
                self.lb = np.min(self.chrom, axis=0, keepdims=True)
            if i == 100:
                self.eps /= 2
            elif i == 500:
                self.eps /= 2
            self.new_chrom = self.chrom.copy()
            self.cross()
            self.mutation()
            self.new_chrom = np.maximum(np.minimum(self.new_chrom, self.ub), self.lb)
            self.choose()
            best_fit.append(np.max(self.fitness))
            mean_fit.append(np.mean(self.fitness))
            # print(np.mean(self.fitness), np.max(self.fitness))
            if i % self.save_freq == 0:
                saved_idx.append(i)
                saved_chrom.append(self.decode())
        saved_idx.append(self.iterations)
        saved_chrom.append(self.decode())
        return np.array(mean_fit), np.array(best_fit), saved_idx, np.array(saved_chrom)


def f(x):
    return np.dot(x, x)


def test_real():
    n_var = 3
    ea = EAReal(n_var, np.ones(n_var) * (-10), np.ones(n_var) * 10, f)
    ea.main()
    # print('-------------------------------')
    # ea = EAReal(n_var, ea.lb, ea.ub, f)
    # ea.main()
    # print('-------------------------------')
    # ea = EAReal(n_var, ea.lb, ea.ub, f)
    # ea.main()
    print(ea.chrom)


def test_bin():
    n_var = 3
    ea = EABi(n_var, np.ones(n_var) * (-10), np.ones(n_var) * 10, f)
    ea.main()
    print(ea.decode())


if __name__ == '__main__':
    test_real()

from ea import *
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def Ball(x):
    # min x^2 + y^2
    # [-5.12, 5.12]
    return x[0] * x[0] + x[1] * x[1]


def DeJong(x):
    # min 100(y - x^2)^2 + (x - 1)^2
    # [-2.048, 2.048]
    return 100 * (x[1] - x[0] ** 2) ** 2 + (x[0] - 1) * (x[0] - 1)


def HimmelBaut(x):
    # min
    # [-6, 6]
    x, y = x[0], x[1]
    return (x * x + y - 11) * (x * x + y - 11) + (x + y * y - 7) * (x + y * y - 7)


def SixHumpCamelBack(x):
    # min
    # [-5, 5]
    # C = 3
    x, y = x[0], x[1]
    return 4 * x * x - 2.1 * x ** 4 + 1 / 3 * x ** 6 + x * y - 4 * y * y + 4 * y ** 4


def Shubert(x):
    # min
    # [-10, 10]
    # C = 190
    x, y = x[0], x[1]
    return (sum([i * np.cos((i + 1) * x + i) for i in range(1, 6)])) * (
        sum([i * np.cos((i + 1) * y + i) for i in range(1, 6)]))


def draw(fun, lb, ub, idx, chrom, title, fig):
    for i in range(5):
        if 'Real' in title:
            ax = fig.add_subplot(6, 5, i + 16)
        else:
            ax = fig.add_subplot(6, 5, i + 21)
        eps = (ub - lb) / 50
        x = np.arange(lb - eps * 5, ub + eps * 5, eps)
        y = np.arange(lb - eps * 5, ub + eps * 5, eps)
        xx, yy = np.meshgrid(x, y)
        zz = fun(np.array([xx, yy]))
        ax.contourf(xx, yy, zz, cmap='rainbow', levels=100)
        x = chrom[i, :, 0]
        y = chrom[i, :, 1]
        ax.scatter(x, y, s=2, marker='x', c='black')
        ax.set_title(title + ' in ' + str(idx[i]), fontsize=8)
        ax.tick_params(labelsize=6)


def wrapper(x):
    return 1 / (1 + x)


def main():
    fun = Ball
    n_var = 2
    lb = np.ones(n_var) * (-5.12)
    ub = np.ones(n_var) * (5.12)

    ea = EAReal(n_var=n_var, lb=lb, ub=ub, ff=fun, wf=wrapper)
    real_mean, real_best, real_id, real_chrom = ea.main()
    real_res = 'EAReal: ' + str(ea.decode()[0]) + ' ' + str(real_best[-1])
    print(real_res)

    ea = EABi(n_var=n_var, lb=lb, ub=ub, ff=fun, wf=wrapper)
    bi_mean, bi_best, bi_id, bi_chrom = ea.main()
    bi_res = 'EABi: ' + str(ea.decode()[0]) + ' ' + str(bi_best[-1])
    print(bi_res)

    fig = plt.figure(1, figsize=(10, 12))

    ax = fig.add_subplot(2, 1, 1, projection='3d')
    eps = (ub[0] - lb[0]) / 100
    x = np.arange(lb[0] - eps * 10, ub[0] + eps * 10, eps)
    y = np.arange(lb[0] - eps * 10, ub[0] + eps * 10, eps)
    xx, yy = np.meshgrid(x, y)
    zz = fun(np.array([xx, yy]))
    ax.plot_surface(xx, yy, zz, cmap='rainbow')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('g(x, y)')
    ax.set_title('3D plot')

    draw(fun, lb[0], ub[0], real_id, real_chrom, 'EAReal', fig)
    draw(fun, lb[0], ub[0], bi_id, bi_chrom, 'EABi', fig)
    ax1 = fig.add_subplot(616)
    idx = np.arange(0, 2000, 10, dtype=int)
    plt1 = ax1.plot(idx, real_mean[idx], 'r', label='EAReal-mean')
    plt2 = ax1.plot(idx, bi_mean[idx], 'g', label='EABi-mean')
    ax2 = ax1.twinx()
    plt3 = ax2.plot(real_best, '--', color='orange', label='EAReal-best')
    plt4 = ax2.plot(bi_best, '--', color='b', label='EABi-best')
    line = plt1 + plt3 + plt2 + plt4
    ax1.set_xlabel('iterations')
    ax1.set_ylabel('mean fitness')
    ax2.set_ylabel('best fitness')
    plt.legend(line, [l.get_label() for l in line], fontsize=7)
    plt.subplots_adjust(left=None, bottom=None, right=None, top=None,
                        wspace=0.5, hspace=0.4)
    plt.text(0, -1, real_res + 'n' + bi_res, size=12, transform=ax1.transAxes)
    plt.savefig('2' + str(fun.__name__) + '.pdf', bbox_inches='tight')


if __name__ == '__main__':
    main()

from ea import *
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def Ball(x):
    # min x^2 + y^2
    # [-100, 100]
    return np.sum(x * x, axis=0)


def Step(x):
    # min
    # [-10, 10]
    return np.sum(np.floor(x + 0.5) ** 2, axis=0)


def Rastrigin(x):
    # min
    # [-5.12, 5.12]
    A = 10
    return np.sum(A + x * x - A * np.cos(2 * np.pi * x), axis=0)


def RosenBrock(x):
    # min
    # [-2.048, 2.048]
    x0 = x[:-1]
    x1 = x[1:]
    return np.sum((1 - x0) ** 2 + 100 * (x1 - x0 ** 2) ** 2, axis=0)


def Ackley(x):
    # min
    # [-32, 32]
    s1 = np.sum(x * x, axis=0)
    s2 = np.sum(np.cos(2 * np.pi * x), axis=0)
    n = x.shape[0]
    return -20 * np.exp(-0.2 * np.sqrt(s1 / n)) - np.exp(s2 / n) + 20 + np.e


def draw(fun, lb, ub, idx, chrom, title, fig):
    for i in range(5):
        if 'Real' in title:
            ax = fig.add_subplot(6, 5, i + 16)
        else:
            ax = fig.add_subplot(6, 5, i + 21)
        eps = (ub - lb) / 100
        x = np.arange(lb - eps * 10, ub + eps * 10, eps)
        y = np.arange(lb - eps * 10, ub + eps * 10, eps)
        xx, yy = np.meshgrid(x, y)
        zz = fun(np.array([xx, yy]))
        ax.contourf(xx, yy, zz, cmap='rainbow', levels=100)
        x = chrom[i, :, 0]
        y = chrom[i, :, 1]
        ax.scatter(x, y, s=2, marker='x', c='black')
        ax.set_title(title + ' in ' + str(idx[i]), fontsize=8)
        ax.tick_params(labelsize=6)


def wrapper(x):
    return 1 / (1 + x)


def main():
    fun = Ball
    n_var = 10
    lb = np.ones(n_var) * (-100)
    ub = np.ones(n_var) * (100)

    ea = EAReal(n_var=n_var, lb=lb, ub=ub, ff=fun, wf=wrapper)
    real_mean, real_best, real_id, real_chrom = ea.main()
    real_res = 'EAReal: ' + str(ea.decode()[0]) + ' ' + str(real_best[-1])
    print(real_res)

    ea = EABi(n_var=n_var, lb=lb, ub=ub, ff=fun, wf=wrapper)
    bi_mean, bi_best, bi_id, bi_chrom = ea.main()
    bi_res = 'EABi: ' + str(ea.decode()[0]) + ' ' + str(bi_best[-1])
    print(bi_res)

    fig = plt.figure(1, figsize=(10, 12))

    ax = fig.add_subplot(2, 1, 1, projection='3d')
    eps = (ub[0] - lb[0]) / 200
    x = np.arange(lb[0] - eps * 10, ub[0] + eps * 10, eps)
    y = np.arange(lb[0] - eps * 10, ub[0] + eps * 10, eps)
    xx, yy = np.meshgrid(x, y)
    zz = fun(np.array([xx, yy]))
    ax.plot_surface(xx, yy, zz, cmap='rainbow')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('g(x, y)')
    ax.set_title('3D plot')

    draw(fun, lb[0], ub[0], real_id, real_chrom, 'EAReal', fig)
    draw(fun, lb[0], ub[0], bi_id, bi_chrom, 'EABi', fig)
    ax1 = fig.add_subplot(616)
    idx = np.arange(0, 2000, 10, dtype=int)
    plt1 = ax1.plot(idx, real_mean[idx], 'r', label='EAReal-mean')
    plt2 = ax1.plot(idx, bi_mean[idx], 'g', label='EABi-mean')
    ax2 = ax1.twinx()
    plt3 = ax2.plot(real_best, '--', color='orange', label='EAReal-best')
    plt4 = ax2.plot(bi_best, '--', color='b', label='EABi-best')
    line = plt1 + plt3 + plt2 + plt4
    ax1.set_xlabel('iterations')
    ax1.set_ylabel('mean fitness')
    ax2.set_ylabel('best fitness')
    plt.legend(line, [l.get_label() for l in line], fontsize=7)
    plt.subplots_adjust(left=None, bottom=None, right=None, top=None,
                        wspace=0.5, hspace=0.4)
    plt.text(0, -1, real_res + 'n' + bi_res, size=12, transform=ax1.transAxes)
    plt.savefig('n' + str(fun.__name__) + '.pdf', bbox_inches='tight')


if __name__ == '__main__':
    main()