HTSim/ga_new.py

import numpy
import numpy as np


class GeneticAlgorithm:
    """Genetic algorithm.
    Implementation of the real-valued Genetic algorithm. The mutations are
    normally distributed perturbations, the selection mechanism is a tournament
    selection, and the crossover oepration is the standard linear combination
    taken at a randomly generated cutting point.
    The total number of evaluations are popsize x ngen
    :param function: Function that can be used to evaluate the entire
        population. It needs to take an input of size pop_size x dim and
        return a numpy.array of size pop_size x 1
    :type function: Object
    :param dim: Number of dimensions
    :type dim: int
    :param lb: Lower variable bounds, of length dim
    :type lb: numpy.array
    :param ub: Lower variable bounds, of length dim
    :type ub: numpy.array
    :param int_var: List of indices with the integer valued variables
        (e.g., [0, 1, 5])
    :type int_var: list
    :param pop_size: Population size
    :type pop_size: int
    :param num_gen: Number of generations
    :type num_gen: int
    :param start: Method for generating the initial population
    :type start: string
    :ivar nvariables: Number of variables (dimensions)
    :ivar nindividuals: population size
    :ivar lower_boundary: lower bounds for the optimization problem
    :ivar upper_boundary: upper bounds for the optimization problem
    :ivar integer_variables: List of variables that are integer valued
    :ivar start: Method for generating the initial population
    :ivar sigma: Perturbation radius. Each pertubation is N(0, sigma)
    :ivar p_mutation: Mutation probability (1/dim)
    :ivar tournament_size: Size of the tournament (5)
    :ivar p_cross: Cross-over probability (0.9)
    :ivar ngenerations: Number of generations
    :ivar function: Object that can be used to evaluate the objective function
    """

    def __init__(self, function, dim, lb, ub, int_var=None, pop_size=6, num_gen=300, start="Random"):

        self.nvariables = dim  # column
        self.nindividuals = pop_size + (pop_size % 2)  # Make sure this is even row
        self.lower_boundary = np.array(lb)
        self.upper_boundary = np.array(ub)
        self.integer_variables = []
        if int_var is not None:
            self.integer_variables = np.array(int_var)
        self.start = start
        self.sigma = 0.2
        self.p_mutation = 1.0 / dim
        self.tournament_size = 5
        self.p_cross = 0.9
        self.ngenerations = num_gen
        self.function = function

    def optimize(self):
        """Method used to run the Genetic algorithm
        :return: Returns the best individual and its function value
        :rtype: numpy.array, float
        """
        #  Initialize population
        if isinstance(self.start, np.ndarray):
            if self.start.shape[0] != self.nindividuals or self.start.shape[1] != self.nvariables:
                raise ValueError("Initial population has incorrect size")
            if any(np.min(self.start, axis=0) < self.lower_boundary) or any(
                    np.max(self.start, axis=0) > self.upper_boundary
            ):
                raise ValueError("Initial population is outside the domain", self.lower_boundary, self.upper_boundary,
                                 self.start)
            population = self.start
        elif self.start == "SLHD":
            from pySOT.experimental_design import SymmetricLatinHypercube

            exp_des = SymmetricLatinHypercube(self.nvariables, self.nindividuals)
            population = self.lower_boundary + exp_des.generate_points() * (self.upper_boundary - self.lower_boundary)
        elif self.start == "LHD":
            from pySOT.experimental_design import LatinHypercube

            exp_des = LatinHypercube(self.nvariables, self.nindividuals)
            population = self.lower_boundary + exp_des.generate_points() * (self.upper_boundary - self.lower_boundary)
        elif self.start == "Random":
            population = self.lower_boundary + np.random.rand(self.nindividuals, self.nvariables) * (
                    self.upper_boundary - self.lower_boundary
            )
        else:
            raise ValueError("Unknown argument for initial population")

        new_population = []
        #  Round positions
        if len(self.integer_variables) > 0:  # 对特定列进行操作
            new_population = np.copy(population)
            population[:, self.integer_variables] = np.round(population[:, self.integer_variables])
            for i in self.integer_variables:
                ind = np.where(population[:, i] < self.lower_boundary[i])
                population[ind, i] += 1
                ind = np.where(population[:, i] > self.upper_boundary[i])
                population[ind, i] -= 1

        for pop in population:
            for i in range(len(pop) // 2):
                if pop[i] >= pop[i + len(pop) // 2]:
                    pop[i], pop[i + len(pop) // 2] = pop[i + len(pop) // 2], pop[i]

        #  Evaluate all individuals
        # function_values = self.function(population) we cannot compute in this way to ensure x is one-dim in policy
        n_row, n_dim = population.shape
        function_values = []
        for r in range(n_row):
            function_values.append(self.function(population[r, :]))
        function_values = np.array(function_values)

        if len(function_values.shape) == 2:
            function_values = np.squeeze(np.asarray(function_values))

        # Save the best individual
        ind = np.argmin(function_values)
        best_individual = np.copy(population[ind, :])  # 找到最优个体
        best_value = function_values[ind]

        if len(self.integer_variables) > 0:
            population = new_population

        # Main loop
        for _ in range(self.ngenerations):
            print('------------------------------')
            print("当前为第{}代".format(_))
            print("最优s为:{}".format(best_individual[0:115]))
            print("最优S为:{}".format(best_individual[115:230]))
            print("最优值为:{}".format(best_value))
            print("------------------------------")
            # Do tournament selection to select the parents
            competitors = np.random.randint(0, self.nindividuals, (self.nindividuals, self.tournament_size))
            ind = np.argmin(function_values[competitors], axis=1)
            winner_indices = np.zeros(self.nindividuals, dtype=int)
            for i in range(self.tournament_size):  # This loop is short
                winner_indices[np.where(ind == i)] = competitors[np.where(ind == i), i]

            parent1 = population[winner_indices[0: self.nindividuals // 2], :]
            parent2 = population[winner_indices[self.nindividuals // 2: self.nindividuals], :]

            # Averaging Crossover
            cross = np.where(np.random.rand(self.nindividuals // 2) < self.p_cross)[0]
            nn = len(cross)  # Number of crossovers
            alpha = np.random.rand(nn, 1)

            # Create the new chromosomes
            parent1_new = np.multiply(alpha, parent1[cross, :]) + np.multiply(1 - alpha, parent2[cross, :])
            parent2_new = np.multiply(alpha, parent2[cross, :]) + np.multiply(1 - alpha, parent1[cross, :])
            parent1[cross, :] = parent1_new
            parent2[cross, :] = parent2_new
            population = np.concatenate((parent1, parent2))

            # Apply mutation
            scale_factors = self.sigma * (self.upper_boundary - self.lower_boundary)  # Scale
            perturbation = np.random.randn(self.nindividuals, self.nvariables)  # Generate perturbations
            perturbation = np.multiply(perturbation, scale_factors)  # Scale accordingly
            perturbation = np.multiply(
                perturbation, (np.random.rand(self.nindividuals, self.nvariables) < self.p_mutation)
            )

            perturbation = round_vars(perturbation, self.integer_variables, self.lower_boundary, self.upper_boundary)
            population = round_vars(population, self.integer_variables, self.lower_boundary, self.upper_boundary)

            population += perturbation  # Add perturbation
            population = np.maximum(np.reshape(self.lower_boundary, (1, self.nvariables)), population)
            population = np.minimum(np.reshape(self.upper_boundary, (1, self.nvariables)), population)

            # Round chromosomes
            new_population = []
            if len(self.integer_variables) > 0:
                new_population = np.copy(population)
                population = round_vars(population, self.integer_variables, self.lower_boundary, self.upper_boundary)

            # Keep the best individual
            population[0, :] = best_individual

            for pop in population:
                for i in range(len(pop) // 2):
                    if pop[i] >= pop[i + len(pop) // 2]:
                        pop[i], pop[i + len(pop) // 2] = pop[i + len(pop) // 2], pop[i]

            #  Evaluate all individuals
            # function_values = self.function(population) we cannot compute in this way to ensure x is one-dim in policy
            n_row, n_dim = population.shape
            function_values = []
            for r in range(n_row):
                function_values.append(self.function(population[r, :]))
            function_values = np.array(function_values)

            if len(function_values.shape) == 2:
                function_values = np.squeeze(np.asarray(function_values))

            # Save the best individual
            ind = np.argmin(function_values)
            best_individual = np.copy(population[ind, :])
            best_value = function_values[ind]

            # Use the positions that are not rounded
            if len(self.integer_variables) > 0:
                population = new_population

        # return best_individual, best_value


def round_vars(x: np.ndarray, int_var, lb, ub):
    """Round integer variables to closest integer in the domain.
    :param x: Set of points, of size npts x dim
    :type x: numpy.array
    :param int_var: Set of indices of integer variables
    :type int_var: numpy.array
    :param lb: Lower bounds, of size 1 x dim
    :type lb: numpy.array
    :param ub: Upper bounds, of size 1 x dim
    :type ub: numpy.array
    :return: The set of points with the integer variables
        rounded to the closest integer in the domain
    :rtype: numpy.array
    """
    # Make sure we don't violate the bound constraints
    for i in int_var:
        ind = np.where(x[:, i] < lb[i])
        x[ind, i] = lb[i]
        ind = np.where(x[:, i] > ub[i])
        x[ind, i] = ub[i]

    if len(int_var) > 0:
        # Round the original ranged integer variables
        x[:, int_var] = np.round(x[:, int_var])
        x = x.astype(numpy.int32, copy=True)
    else:
        x = x.astype(numpy.float64, copy=True)
    return x