# !/usr/bin/env python # -*- coding:utf-8 -*- # Author: wsw # 簡單實現SGA算法 import numpy as np from scipy.optimize import fsolve, basinhopping import random import timeit   # 根據解的精度確定染色體(chromosome)的長度 # 需要根據決策變量的上下邊界來確定 def getEncodedLength(delta=0.0001, boundarylist=[]):  # 每個變量的編碼長度  lengths = []  for i in boundarylist:   lower = i[0]   upper = i[1]   # lamnda 代表匿名函數f(x)=0,50代表搜索的初始解   res = fsolve(lambda x: ((upper - lower) * 1 / delta) - 2 ** x - 1, 50)   length = int(np.floor(res[0]))   lengths.append(length)  return lengths  pass   # 隨機生成初始編碼種群 def getIntialPopulation(encodelength, populationSize):  # 隨機化初始種群為0  chromosomes = np.zeros((populationSize, sum(encodelength)), dtype=np.uint8)  for i in range(populationSize):   chromosomes[i, :] = np.random.randint(0, 2, sum(encodelength))  # print('chromosomes shape:', chromosomes.shape)  return chromosomes   # 染色體解碼得到表現型的解 def decodedChromosome(encodelength, chromosomes, boundarylist, delta=0.0001):  populations = chromosomes.shape[0]  variables = len(encodelength)  decodedvalues = np.zeros((populations, variables))  for k, chromosome in enumerate(chromosomes):   chromosome = chromosome.tolist()   start = 0   for index, length in enumerate(encodelength):    # 將一個染色體進行拆分，得到染色體片段    power = length - 1    # 解碼得到的10進制數字    demical = 0    for i in range(start, length + start):     demical += chromosome[i] * (2 ** power)     power -= 1    lower = boundarylist[index][0]    upper = boundarylist[index][1]    decodedvalue = lower + demical * (upper - lower) / (2 ** length - 1)    decodedvalues[k, index] = decodedvalue    # 開始去下一段染色體的編碼    start = length  return decodedvalues   # 得到個體的適應度值及每個個體被選擇的累積概率 def getFitnessValue(func, chromosomesdecoded):  # 得到種群規模和決策變量的個數  population, nums = chromosomesdecoded.shape  # 初始化種群的適應度值為0  fitnessvalues = np.zeros((population, 1))  # 計算適應度值  for i in range(population):   fitnessvalues[i, 0] = func(chromosomesdecoded[i, :])  # 計算每個染色體被選擇的概率  probability = fitnessvalues / np.sum(fitnessvalues)  # 得到每個染色體被選中的累積概率  cum_probability = np.cumsum(probability)  return fitnessvalues, cum_probability   # 新種群選擇 def selectNewPopulation(chromosomes, cum_probability):  m, n = chromosomes.shape  newpopulation = np.zeros((m, n), dtype=np.uint8)  # 隨機產生M個概率值  randoms = np.random.rand(m)  for i, randoma in enumerate(randoms):   logical = cum_probability >= randoma   index = np.where(logical == 1)   # index是tuple,tuple中元素是ndarray   newpopulation[i, :] = chromosomes[index[0][0], :]  return newpopulation  pass   # 新種群交叉 def crossover(population, Pc=0.8):  """  :param population: 新種群  :param Pc: 交叉概率默認是0.8  :return: 交叉后得到的新種群  """  # 根據交叉概率計算需要進行交叉的個體個數  m, n = population.shape  numbers = np.uint8(m * Pc)  # 確保進行交叉的染色體個數是偶數個  if numbers % 2 != 0:   numbers += 1  # 交叉后得到的新種群  updatepopulation = np.zeros((m, n), dtype=np.uint8)  # 產生隨機索引  index = random.sample(range(m), numbers)  # 不進行交叉的染色體進行復制  for i in range(m):   if not index.__contains__(i):    updatepopulation[i, :] = population[i, :]  # crossover  while len(index) > 0:   a = index.pop()   b = index.pop()   # 隨機產生一個交叉點   crossoverPoint = random.sample(range(1, n), 1)   crossoverPoint = crossoverPoint[0]   # one-single-point crossover   updatepopulation[a, 0:crossoverPoint] = population[a, 0:crossoverPoint]   updatepopulation[a, crossoverPoint:] = population[b, crossoverPoint:]   updatepopulation[b, 0:crossoverPoint] = population[b, 0:crossoverPoint]   updatepopulation[b, crossoverPoint:] = population[a, crossoverPoint:]  return updatepopulation  pass   # 染色體變異 def mutation(population, Pm=0.01):  """   :param population: 經交叉后得到的種群  :param Pm: 變異概率默認是0.01  :return: 經變異操作后的新種群  """  updatepopulation = np.copy(population)  m, n = population.shape  # 計算需要變異的基因個數  gene_num = np.uint8(m * n * Pm)  # 將所有的基因按照序號進行10進制編碼，則共有m*n個基因  # 隨機抽取gene_num個基因進行基本位變異  mutationGeneIndex = random.sample(range(0, m * n), gene_num)  # 確定每個將要變異的基因在整個染色體中的基因座(即基因的具體位置)  for gene in mutationGeneIndex:   # 確定變異基因位于第幾個染色體   chromosomeIndex = gene // n   # 確定變異基因位于當前染色體的第幾個基因位   geneIndex = gene % n   # mutation   if updatepopulation[chromosomeIndex, geneIndex] == 0:    updatepopulation[chromosomeIndex, geneIndex] = 1   else:    updatepopulation[chromosomeIndex, geneIndex] = 0  return updatepopulation  pass   # 定義適應度函數 def fitnessFunction():  return lambda x: 21.5 + x[0] * np.sin(4 * np.pi * x[0]) + x[1] * np.sin(20 * np.pi * x[1])  pass   def main(max_iter=500):  # 每次迭代得到的最優解  optimalSolutions = []  optimalValues = []  # 決策變量的取值范圍  decisionVariables = [[-3.0, 12.1], [4.1, 5.8]]  # 得到染色體編碼長度  lengthEncode = getEncodedLength(boundarylist=decisionVariables)  for iteration in range(max_iter):   # 得到初始種群編碼   chromosomesEncoded = getIntialPopulation(lengthEncode, 10)   # 種群解碼   decoded = decodedChromosome(lengthEncode, chromosomesEncoded, decisionVariables)   # 得到個體適應度值和個體的累積概率   evalvalues, cum_proba = getFitnessValue(fitnessFunction(), decoded)   # 選擇新的種群   newpopulations = selectNewPopulation(chromosomesEncoded, cum_proba)   # 進行交叉操作   crossoverpopulation = crossover(newpopulations)   # mutation   mutationpopulation = mutation(crossoverpopulation)   # 將變異后的種群解碼，得到每輪迭代最終的種群   final_decoded = decodedChromosome(lengthEncode, mutationpopulation, decisionVariables)   # 適應度評價   fitnessvalues, cum_individual_proba = getFitnessValue(fitnessFunction(), final_decoded)   # 搜索每次迭代的最優解，以及最優解對應的目標函數的取值   optimalValues.append(np.max(list(fitnessvalues)))   index = np.where(fitnessvalues == max(list(fitnessvalues)))   optimalSolutions.append(final_decoded[index[0][0], :])  # 搜索最優解  optimalValue = np.max(optimalValues)  optimalIndex = np.where(optimalValues == optimalValue)  optimalSolution = optimalSolutions[optimalIndex[0][0]]  return optimalSolution, optimalValue   solution, value = main() print('最優解: x1, x2') print(solution[0], solution[1]) print('最優目標函數值:', value) # 測量運行時間 elapsedtime = timeit.timeit(stmt=main, number=1) print('Searching Time Elapsed:(S)', elapsedtime)