请修改这段代码，使得运行过程中指数函数的指数部分只有常数 # -*- coding: UTF-8 -*- import os import pandas as pd import math import random import operator from deap import creator, base, tools, gp, algorithms import numpy from fontTools.misc.py23 import xrange from matplotlib import pyplot as plt import pygraphviz as pgv import datetime import sys path = os.getcwd() curr_time = datetime.datetime.now() time_str = datetime.datetime.strftime(curr_time, '%Y%m%d%H%M') def text_save(filename, data): file = open(filename, 'a') for i in range(len(data)): s = str(data[i]).replace('[', '').replace(']', '').replace('\n', '') s = s.replace("'", '').replace(',', '') + '\n' file.write(s) file.close() print("Saving the file succeeded") return def ProDiv(left, right): if right == 0: return 1000000 else: return left / right def ProSqrt(x): if x <= 0: return 1000000 else: return x ** 0.5 def Prolog(x): if x <= 0: return 1000000 else: return math.log(x) def Proexp(x): if x >= 100: return 1000000 else: return math.exp(x) def Prosquare(x): if x <= 0: return 1000000 else: if x ** 2 > 1e30: return 1000000 else: return x ** 2 def Propow(down, up): """ 安全的指数计算函数，指数必须是纯常数且范围在(-10,10) 参数: down: 底数（可以包含变量） up: 指数（必须是纯常数，范围在-10到10之间） 返回: 计算结果，错误时返回1000000 """ # 确保指数是有效常数 # 1. 检查是否为数字类型 if not isinstance(up, (int, float)): return 1000000 # 2. 检查是否为特殊值（NaN或无穷大） if math.isnan(up) or math.isinf(up): return 1000000 # 3. 检查指数范围（-10 < up < 10） if up <= -10 or up >= 10: return 1000000 # 处理底数为零的情况 if down == 0: return 0.0 if up > 0 else 1000000 # 0的正指数=0，负指数未定义 # 处理负底数 if down < 0: return 1000000 # 使用对数预测结果大小，防止溢出 try: # 计算结果的绝对值的对数 log_value = abs(up) * math.log10(abs(down)) # 检查结果是否过大(>10^50)或过小(<10^-50) if log_value > 50 or log_value < -50: return 1000000 except (ValueError, OverflowError): return 1000000 # 正常计算并捕获可能的溢出错误 try: result = down ** up # 额外检查是否得到无穷大 if math.isinf(result): return 1000000 return result except (OverflowError, ValueError): return 1000000 # 更可靠的常数节点检查 def is_constant_node(node): """检查节点是否是真正的常数节点""" # 常数节点必须有 value 属性 if not hasattr(node, 'value'): return False # 确保值是可接受的数字类型 if not isinstance(node.value, (int, float)): return False # 确保值在合理范围内 if abs(node.value) > 1e10: # 防止过大或过小的常数 return False return True # 更严格的幂函数检查 def check_pow_arguments(individual): """确保所有 Propow 节点的指数参数都是常数""" # 首先将个体转换为字符串表达式 expr_str = str(individual) # 查找所有 Propow 函数调用 pow_calls = [] start_idx = 0 while True: idx = expr_str.find("Propow(", start_idx) if idx == -1: break # 找到匹配的右括号 depth = 1 end_idx = idx + len("Propow(") while end_idx < len(expr_str) and depth > 0: if expr_str[end_idx] == '(': depth += 1 elif expr_str[end_idx] == ')': depth -= 1 end_idx += 1 # 提取参数部分 args_str = expr_str[idx + len("Propow("):end_idx - 1].strip() args = [arg.strip() for arg in args_str.split(',') if arg.strip()] # 确保有两个参数 if len(args) == 2: pow_calls.append((args[0], args[1])) start_idx = end_idx # 检查每个幂调用的指数参数 for base, exponent in pow_calls: # 指数必须是数字（常数） try: float(exponent) except ValueError: # 如果无法转换为浮点数，说明是变量 print(f"发现非法幂函数表达式: Propow({base}, {exponent})") return False return True # 在初始种群生成时强制验证 def create_valid_individual(): """创建并验证个体，直到找到有效的个体""" max_attempts = 100 attempts = 0 while attempts < max_attempts: ind = tool.individual() if check_pow_arguments(ind): return ind attempts += 1 # 如果多次尝试后仍无法生成合法个体，创建一个安全的默认个体 # 例如: Propow(1.0, 1.0) default_expr = gp.PrimitiveTree.from_string( "Propow(1.0, 1.0)", pset ) return creator.Individual(default_expr) # 修改种群创建方法 def create_valid_population(n): """创建并验证种群""" pop = [] for _ in range(n): pop.append(create_valid_individual()) return pop # 增强的交叉操作 def custom_cxOnePoint(ind1, ind2): max_attempts = 50 attempts = 0 while attempts < max_attempts: new_ind1, new_ind2 = gp.cxOnePoint(ind1, ind2) if check_pow_arguments(new_ind1) and check_pow_arguments(new_ind2): return new_ind1, new_ind2 attempts += 1 # 如果多次尝试后仍无法生成合法个体，返回原始个体 return ind1, ind2 # 增强的变异操作 def custom_mutate(individual, expr, pset): max_attempts = 50 attempts = 0 while attempts < max_attempts: new_individual, = gp.mutUniform(individual, expr, pset) if check_pow_arguments(new_individual): return new_individual, attempts += 1 # 如果多次尝试后仍无法生成合法个体，返回原始个体 return individual, def Psin(l, r): return math.sin(l * r) def Pcos(l, r): return math.cos(l * r) def generate_random_float(min_value, max_value, decimal_places): scale = 10 ** decimal_places random_float = random.uniform(min_value, max_value) rounded_float = round(random_float, decimal_places) scaled_float = int(rounded_float * scale) / scale return scaled_float min_val = 0.0 max_val = 1.0 decimal_places = 4 # 定义 primitive set pset = gp.PrimitiveSet('main', 4) # 定义变量个数 pset.renameArguments(ARG0='Fr') # 定义变量名称 pset.renameArguments(ARG1='We') pset.renameArguments(ARG2='Re') pset.renameArguments(ARG3='Bo') # 只保留乘除和幂运算 pset.addPrimitive(operator.mul, 2, name='mul') pset.addPrimitive(ProDiv, 2, name='ProDiv') pset.addPrimitive(Propow, 2, name='Propow') # 添加常数生成器 pset.addEphemeralConstant('c1', lambda: generate_random_float(-2, 2, 4)) pset.addEphemeralConstant('c2', lambda: generate_random_float(-1, 1, 4)) pset.addEphemeralConstant('c3', lambda: generate_random_float(0, 1, 4)) pset.addEphemeralConstant('c4', lambda: generate_random_float(1, 2, 4)) # 创建fitness类、individual类 creator.create('FitnessMin', base.Fitness, weights=(-1.0,)) creator.create('Individual', gp.PrimitiveTree, fitness=creator.FitnessMin) # 定义个体的生成方法、种群的生成方法 tool = base.Toolbox() tool.register('expr', gp.genHalfAndHalf, pset=pset, min_=1, max_=3) tool.register('individual', tools.initIterate, creator.Individual, tool.expr) tool.register('population', tools.initRepeat, list, tool.individual) tool.register('compile', gp.compile, pset=pset) tool.register('expr_mut', gp.genFull, pset=pset, min_=0, max_=2) # 生成一个subtree # 修改变异操作，确保变异后的个体符合要求 def custom_mutate(individual, expr, pset): new_individual, = gp.mutUniform(individual, expr, pset) while not check_pow_arguments(new_individual): new_individual, = gp.mutUniform(individual, expr, pset) return new_individual, # 修改交叉操作，确保交叉后的个体符合要求 def custom_cxOnePoint(ind1, ind2): new_ind1, new_ind2 = gp.cxOnePoint(ind1, ind2) while not check_pow_arguments(new_ind1) or not check_pow_arguments(new_ind2): new_ind1, new_ind2 = gp.cxOnePoint(ind1, ind2) return new_ind1, new_ind2 # 定义评价函数 def fit_evaluation(individual, input1, input2, input3, input4, input5, output): func = tool.compile(expr=individual) # pset上面已经给过了 sqerrors = ((abs((1 + func(Fr, We, Re, Bo)) * dPa - dPe) / dPe * 100) for Fr, We, Re, Bo, dPa, dPe in zip(input1, input2, input3, input4, input5, output)) return math.fsum(sqerrors) / len(output), # 必须返回一个tuple path = os.getcwd() filename = path + '/' + 'data_horizontal_mg_4g' + '.xlsx' book = pd.read_excel(filename) data = book.values # 检查文件是否存在 if os.path.exists(filename): # 确保已经导入了os模块 try: # 读取Excel文件 book = pd.read_excel(filename) data = book.values print(f"成功读取文件: {filename}") print(f"数据形状: {data.shape}") except Exception as e: print(f"读取文件时发生错误: {e}") else: print(f"错误：文件不存在 - {filename}") # 可以在这里添加创建文件或其他处理逻辑 # x = data[:, 37] # X = (data[:, 31]/data[:, 32])**0.5 Fr = data[:, 28] # (data[:, 1]+2.2) We = data[:, 32] Re = data[:, 24] Bo = data[:, 36] # N=(WeG*Fr**2)**0.25 # Pred = data[:, 19] # f = data[:, 35]/data[:, 36] dPl = data[:, 61] dPe = data[:, 65] Y = (data[:, 67] / data[:, 66]) ** 0.5 dPa = data[:, 66] # /1000000000 # 定义evaluate、select、mate、mutate（这几个名字必须这样取，否则出错） tool.register('evaluate', fit_evaluation, input1=Fr, input2=We, input3=Re, input4=Bo, input5=dPa, output=dPe) tool.register('select', tools.selTournament, tournsize=3) # 注册修改后的交叉操作 tool.register('mate', custom_cxOnePoint) # 注册修改后的变异操作 tool.register('mutate', custom_mutate, expr=tool.expr_mut, pset=pset) # 限制树深度 tool.decorate("mate", gp.staticLimit(key=operator.attrgetter("height"), max_value=10)) tool.decorate("mutate", gp.staticLimit(key=operator.attrgetter("height"), max_value=10)) def drawBestInd(expr): # 输入参数是一个表达式（expr/tree），这个例子中就是individual nodes, edges, labels = gp.graph(expr) g = pgv.AGraph() g.add_nodes_from(nodes) g.add_edges_from(edges) g.layout(prog="dot") for i in nodes: n = g.get_node(i) n.attr["label"] = labels[i] # curr_time = datetime.datetime.now() # time_str = datetime.datetime.strftime(curr_time, '%Y%m%d%H%M') g.draw(time_str + "tree.pdf") def varOr(population, toolbox, lambda_, cxpb, mutpb): """Part of an evolutionary algorithm applying only the variation part (crossover, mutation **or** reproduction). The modified individuals have their fitness invalidated. The individuals are cloned so returned population is independent of the input population. :param population: A list of individuals to vary. :param toolbox: A :class:`~deap.base.Toolbox` that contains the evolution operators. :param lambda\_: The number of children to produce :param cxpb: The probability of mating two individuals. :param mutpb: The probability of mutating an individual. :returns: The final population. The variation goes as follow. On each of the *lambda_* iteration, it selects one of the three operations; crossover, mutation or reproduction. In the case of a crossover, two individuals are selected at random from the parental population :math:`P_\mathrm{p}`, those individuals are cloned using the :meth:`toolbox.clone` method and then mated using the :meth:`toolbox.mate` method. Only the first child is appended to the offspring population :math:`P_\mathrm{o}`, the second child is discarded. In the case of a mutation, one individual is selected at random from :math:`P_\mathrm{p}`, it is cloned and then mutated using using the :meth:`toolbox.mutate` method. The resulting mutant is appended to :math:`P_\mathrm{o}`. In the case of a reproduction, one individual is selected at random from :math:`P_\mathrm{p}`, cloned and appended to :math:`P_\mathrm{o}`. This variation is named *Or* because an offspring will never result from both operations crossover and mutation. The sum of both probabilities shall be in :math:`[0, 1]`, the reproduction probability is 1 - *cxpb* - *mutpb*. """ assert (cxpb + mutpb) <= 1.0, ( "The sum of the crossover and mutation probabilities must be smaller " "or equal to 1.0.") offspring = [] for _ in xrange(lambda_): op_choice = random.random() if op_choice < cxpb: # Apply crossover ind1, ind2 = map(toolbox.clone, random.sample(population, 2)) ind1, ind2 = toolbox.mate(ind1, ind2) del ind1.fitness.values offspring.append(ind1) elif op_choice < cxpb + mutpb: # Apply mutation ind = toolbox.clone(random.choice(population)) ind, = toolbox.mutate(ind) del ind.fitness.values offspring.append(ind) else: # Apply reproduction offspring.append(random.choice(population)) return offspring def eaMuPlusLambda(population, toolbox, mu, lambda_, cxpb, mutpb, ngen, stats=None, halloffame=None, verbose=__debug__): """This is the :math:`(\mu + \lambda)` evolutionary algorithm. :param population: A list of individuals. :param toolbox: A :class:`~deap.base.Toolbox` that contains the evolution operators. :param mu: The number of individuals to select for the next generation. :param lambda\_: The number of children to produce at each generation. :param cxpb: The probability that an offspring is produced by crossover. :param mutpb: The probability that an offspring is produced by mutation. :param ngen: The number of generation. :param stats: A :class:`~deap.tools.Statistics` object that is updated inplace, optional. :param halloffame: A :class:`~deap.tools.HallOfFame` object that will contain the best individuals, optional. :param verbose: Whether or not to log the statistics. :returns: The final population :returns: A class:`~deap.tools.Logbook` with the statistics of the evolution. The algorithm takes in a population and evolves it in place using the :func:`varOr` function. It returns the optimized population and a :class:`~deap.tools.Logbook` with the statistics of the evolution. The logbook will contain the generation number, the number of evaluations for each generation and the statistics if a :class:`~deap.tools.Statistics` is given as argument. The *cxpb* and *mutpb* arguments are passed to the :func:`varOr` function. The pseudocode goes as follow :: evaluate(population) for g in range(ngen): offspring = varOr(population, toolbox, lambda_, cxpb, mutpb) evaluate(offspring) population = select(population + offspring, mu) First, the individuals having an invalid fitness are evaluated. Second, the evolutionary loop begins by producing *lambda_* offspring from the population, the offspring are generated by the :func:`varOr` function. The offspring are then evaluated and the next generation population is selected from both the offspring **and** the population. Finally, when *ngen* generations are done, the algorithm returns a tuple with the final population and a :class:`~deap.tools.Logbook` of the evolution. This function expects :meth:`toolbox.mate`, :meth:`toolbox.mutate`, :meth:`toolbox.select` and :meth:`toolbox.evaluate` aliases to be registered in the toolbox. This algorithm uses the :func:`varOr` variation. """ logbook = tools.Logbook() logbook.header = ['gen', 'nevals'] + (stats.fields if stats else []) # Evaluate the individuals with an invalid fitness invalid_ind = [ind for ind in population if not ind.fitness.valid] fitnesses = toolbox.map(toolbox.evaluate, invalid_ind) for ind, fit in zip(invalid_ind, fitnesses): ind.fitness.values = fit if halloffame is not None: halloffame.update(population) record = stats.compile(population) if stats is not None else {} logbook.record(gen=0, nevals=len(invalid_ind), **record) if verbose: print(logbook.stream) # Begin the generational process for gen in range(1, ngen + 1): # Vary the population offspring = varOr(population, toolbox, lambda_, cxpb, mutpb) # Evaluate the individuals with an invalid fitness invalid_ind = [ind for ind in offspring if not ind.fitness.valid] fitnesses = toolbox.map(toolbox.evaluate, invalid_ind) for ind, fit in zip(invalid_ind, fitnesses): ind.fitness.values = fit # Update the hall of fame with the generated individuals if halloffame is not None: halloffame.update(offspring) # Select the next generation population population[:] = toolbox.select(population + offspring, mu) best_ind = halloffame.items[0] # Update the statistics with the new population record = stats.compile(population) if stats is not None else {} logbook.record(gen=gen, nevals=len(invalid_ind), **record) temp = logbook.stream if verbose: print(temp) print(best_ind) f = open(time_str + "outputfile.dat", "a") f.write(str(temp) + '\n') f.write(str(best_ind) + '\n') f.close() # 关闭文件 if halloffame.items[0].fitness.values[0] < 20: print("绝对百分比误差小于20%，完成计算！") print(best_ind) return population, logbook return population, logbook def varAnd(population, toolbox, cxpb, mutpb): r"""Part of an evolutionary algorithm applying only the variation part (crossover **and** mutation). The modified individuals have their fitness invalidated. The individuals are cloned so returned population is independent of the input population. :param population: A list of individuals to vary. :param toolbox: A :class:`~deap.base.Toolbox` that contains the evolution operators. :param cxpb: The probability of mating two individuals. :param mutpb: The probability of mutating an individual. :returns: A list of varied individuals that are independent of their parents. The variation goes as follow. First, the parental population :math:`P_\mathrm{p}` is duplicated using the :meth:`toolbox.clone` method and the result is put into the offspring population :math:`P_\mathrm{o}`. A first loop over :math:`P_\mathrm{o}` is executed to mate pairs of consecutive individuals. According to the crossover probability *cxpb*, the individuals :math:`\mathbf{x}_i` and :math:`\mathbf{x}_{i+1}` are mated using the :meth:`toolbox.mate` method. The resulting children :math:`\mathbf{y}_i` and :math:`\mathbf{y}_{i+1}` replace their respective parents in :math:`P_\mathrm{o}`. A second loop over the resulting :math:`P_\mathrm{o}` is executed to mutate every individual with a probability *mutpb*. When an individual is mutated it replaces its not mutated version in :math:`P_\mathrm{o}`. The resulting :math:`P_\mathrm{o}` is returned. This variation is named *And* because of its propensity to apply both crossover and mutation on the individuals. Note that both operators are not applied systematically, the resulting individuals can be generated from crossover only, mutation only, crossover and mutation, and reproduction according to the given probabilities. Both probabilities should be in :math:`[0, 1]`. """ offspring = [toolbox.clone(ind) for ind in population] # Apply crossover and mutation on the offspring for i in range(1, len(offspring), 2): if random.random() < cxpb: offspring[i - 1], offspring[i] = toolbox.mate(offspring[i - 1], offspring[i]) del offspring[i - 1].fitness.values, offspring[i].fitness.values for i in range(len(offspring)): if random.random() < mutpb: offspring[i], = toolbox.mutate(offspring[i]) del offspring[i].fitness.values return offspring def eaSimple(population, toolbox, cxpb, mutpb, ngen, stats=None, halloffame=None, verbose=__debug__): """This algorithm reproduce the simplest evolutionary algorithm as presented in chapter 7 of [Back2000]_. :param population: A list of individuals. :param toolbox: A :class:`~deap.base.Toolbox` that contains the evolution operators. :param cxpb: The probability of mating two individuals. :param mutpb: The probability of mutating an individual. :param ngen: The number of generation. :param stats: A :class:`~deap.tools.Statistics` object that is updated inplace, optional. :param halloffame: A :class:`~deap.tools.HallOfFame` object that will contain the best individuals, optional. :param verbose: Whether or not to log the statistics. :returns: The final population :returns: A class:`~deap.tools.Logbook` with the statistics of the evolution The algorithm takes in a population and evolves it in place using the :meth:`varAnd` method. It returns the optimized population and a :class:`~deap.tools.Logbook` with the statistics of the evolution. The logbook will contain the generation number, the number of evaluations for each generation and the statistics if a :class:`~deap.tools.Statistics` is given as argument. The *cxpb* and *mutpb* arguments are passed to the :func:`varAnd` function. The pseudocode goes as follow :: evaluate(population) for g in range(ngen): population = select(population, len(population)) offspring = varAnd(population, toolbox, cxpb, mutpb) evaluate(offspring) population = offspring As stated in the pseudocode above, the algorithm goes as follow. First, it evaluates the individuals with an invalid fitness. Second, it enters the generational loop where the selection procedure is applied to entirely replace the parental population. The 1:1 replacement ratio of this algorithm **requires** the selection procedure to be stochastic and to select multiple times the same individual, for example, :func:`~deap.tools.selTournament` and :func:`~deap.tools.selRoulette`. Third, it applies the :func:`varAnd` function to produce the next generation population. Fourth, it evaluates the new individuals and compute the statistics on this population. Finally, when *ngen* generations are done, the algorithm returns a tuple with the final population and a :class:`~deap.tools.Logbook` of the evolution. .. note:: Using a non-stochastic selection method will result in no selection as the operator selects *n* individuals from a pool of *n*. This function expects the :meth:`toolbox.mate`, :meth:`toolbox.mutate`, :meth:`toolbox.select` and :meth:`toolbox.evaluate` aliases to be registered in the toolbox. .. [Back2000] Back, Fogel and Michalewicz, "Evolutionary Computation 1 : Basic Algorithms and Operators", 2000. """ logbook = tools.Logbook() logbook.header = ['gen', 'nevals'] + (stats.fields if stats else []) # Evaluate the individuals with an invalid fitness invalid_ind = [ind for ind in population if not ind.fitness.valid] fitnesses = toolbox.map(toolbox.evaluate, invalid_ind) for ind, fit in zip(invalid_ind, fitnesses): ind.fitness.values = fit if halloffame is not None: halloffame.update(population) record = stats.compile(population) if stats else {} logbook.record(gen=0, nevals=len(invalid_ind), **record) if verbose: print(logbook.stream) # Begin the generational process for gen in range(1, ngen + 1): # Select the next generation individuals offspring = toolbox.select(population, len(population)) # Vary the pool of individuals offspring = varAnd(offspring, toolbox, cxpb, mutpb) # Evaluate the individuals with an invalid fitness invalid_ind = [ind for ind in offspring if not ind.fitness.valid] fitnesses = toolbox.map(toolbox.evaluate, invalid_ind) for ind, fit in zip(invalid_ind, fitnesses): ind.fitness.values = fit # Update the hall of fame with the generated individuals if halloffame is not None: halloffame.update(offspring) # Replace the current population by the offspring population[:] = offspring # Append the current generation statistics to the logbook best_ind = halloffame.items[0] record = stats.compile(population) if stats else {} logbook.record(gen=gen, nevals=len(invalid_ind), **record) temp = logbook.stream if verbose: print(temp) print(best_ind) f = open(time_str + "outputfile.dat", "a") f.write(str(temp) + '\n') f.write(str(best_ind) + '\n') f.close() # 关闭文件 if halloffame.items[0].fitness.values[0] < 20: print("绝对百分比误差小于20%，完成计算！") print(best_ind) return population, logbook return population, logbook def main(): random.seed(318) ### 参数设置 cxpb = 0.5 # 交叉概率 mutpb = 0.15 # 变异概率0 ngen = 1500 # 迭代次数 popSize = 150 # 种群规模 ### 生成初始种群 - 使用新的验证方法 print("创建初始种群...") pop = create_valid_population(popSize) print("初始种群创建完成") # 验证初始种群中的所有个体 for i, ind in enumerate(pop): if not check_pow_arguments(ind): print(f"警告: 初始种群中的个体 {i} 无效!") print(str(ind)) hof = tools.HallOfFame(1) invalid_fit = [ind for ind in pop if not ind.fitness.valid] fitnesses = tool.map(tool.evaluate, invalid_fit) # 求初始种群的每个个体的适应度值，是一个list for ind, fit in zip(invalid_fit, fitnesses): ind.fitness.values = fit # 给每个ind的fitness赋值 hof.update(pop) best_inds = [] # 记录每一代的最优个体 best_ind = hof.items[0] best_inds.append(best_ind) stats = tools.Statistics(lambda ind: ind.fitness.values) stats.register('avg', numpy.mean, axis=0) ### axis = 0表示数组的列求平均；axis = 1 表示数组的行求平均 stats.register('std', numpy.std, axis=0) stats.register('min', numpy.min, axis=0) stats.register('max', numpy.max, axis=0) ### record是一个字典：record{‘avg’：[],'std':[],...},每个value是一个list record = stats.compile(pop) for key in record: record[key] = record[key][0] logbook = tools.Logbook() logbook.record(gen=0, eval=popSize, best_ind=best_ind, **record) logbook.header = 'gen', 'min', 'max', 'avg', 'std', 'eval', 'best_ind' print(logbook) print('--------开始迭代--------') for g in range(ngen): ### select 选出popSize个个体进入下一代 # print('The', g + 1, 'step') offspring = tool.select(pop, len(pop)) # offspring = list(map(tool.clone, offspring)) # offspring = tool.map(tool.clone, offspring) offspring = [tool.clone(ind) for ind in offspring] # Apply crossover and mutation on the offspring for i in range(1, len(offspring), 2): if random.random() < cxpb: offspring[i - 1], offspring[i] = tool.mate(offspring[i - 1], offspring[i]) del offspring[i - 1].fitness.values, offspring[i].fitness.values for i in range(len(offspring)): if random.random() < mutpb: offspring[i], = tool.mutate(offspring[i]) del offspring[i].fitness.values # Evaluate the individuals with an invalid fitness invalid_fit = [ind for ind in offspring if not ind.fitness.valid] fitnesses = tool.map(tool.evaluate, invalid_fit) for ind, fit in zip(invalid_fit, fitnesses): ind.fitness.values = fit hof.update(offspring) # 找到本代的最优个体 pop[:] = offspring # 更新种群 best_ind = hof.items[0] best_inds.append(best_ind) record = stats.compile(pop) # 数据统计 for key in record: record[key] = record[key][0] logbook.record(gen=g + 1, eval=len(invalid_fit), best_ind=best_ind, **record) temp = logbook.stream print(temp) # print logbook.stream ff = open(time_str + "outputfile.dat", "a") ff.write(str(temp) + '\n') # f.write(str(best_ind) + '\n') ff.close() # 关闭文件 if hof.items[0].fitness.values[0] < 20: print("绝对百分比误差小于20%，完成计算！") print(best_ind) print('--------迭代结束-------') break print('--------迭代结束-------') # print(logbook) # pop = tool.population(n=150) # hof = tools.HallOfFame(1) # # stats_fit = tools.Statistics(lambda ind: ind.fitness.values) # stats_size = tools.Statistics(len) # mstats = tools.MultiStatistics(fitness=stats_fit, size=stats_size) # mstats.register("avg", numpy.mean) # mstats.register("std", numpy.std) # mstats.register("min", numpy.min) # mstats.register("max", numpy.max) # pop, log = eaSimple(pop, tool, 0.5, 0.1, 5000, stats=mstats, halloffame=hof, verbose=True) # pop, log = eaMuPlusLambda(pop, tool, 20, 40, 0.5, 0.1, 5000, # stats=mstats, halloffame=hof, verbose=True) # print log # return pop, log, hof drawBestInd(expr=hof.items[0]) # 画出实际函数和优化的函数 # temp = numpy.zeros((len(data))) temp = numpy.zeros((len(data), 5)) func = tool.compile(expr=hof.items[0]) for i in range(len(data)): temp[i, 1] = (func(Fr[i], We[i], Re[i], Bo[i]) + 1) * dPa[ i] # func(FrG, Bo, Rel, ReG, dPa) func(Fr, Bo, f, Pred, x, N) x1 = numpy.arange(0, (max(data[:, 65])), 10) y1 = 1.3 * x1 y2 = 0.7 * x1 # func = tool.compile(expr=hof.items[0]) # z = [func(x) for x in points] temp[:, 0] = data[:, 65] temp[:, 2] = data[:, 66] # temp[:, 3] = data[:, 7] # temp[:, 4] = data[:, 8] filename = time_str + '_data.dat' text_save(filename, temp) # f1 = open(time_str + "outputfile.dat", "a") # f1.write( + '\n') # f1.close ax1 = plt.subplot(211) plt.plot(data[:, 65], temp[:, 1], 'o', markeredgecolor='b', markerfacecolor='none', label='EXP-PRE-ENV') plt.plot(x1, y1, '--', color='k') # , label='+30%') plt.plot(x1, y2, '--', color='k') # , label='-30%') plt.legend(loc='best') plt.xscale('log') plt.yscale('log') ax1 = plt.subplot(212) plt.plot(data[:, 65], data[:, 66], 'o', markeredgecolor='r', markerfacecolor='none', label='EXP-PRE_ANA') plt.plot(x1, y1, '--', color='k') # , label='+30%') plt.plot(x1, y2, '--', color='k') # , label='-30%') plt.legend(loc='best') plt.xscale('log') plt.yscale('log') # ax1 = plt.subplot(223) # plt.plot(data[:, 5], data[:, 7], 'o', markeredgecolor='c', markerfacecolor='none', label='EXP-PRE_MSH') # plt.plot(x1, y1, '--', color='k')#, label='+30%') # plt.plot(x1, y2, '--', color='k')#, label='-30%') # plt.legend(loc='best') # plt.xscale('log') # plt.yscale('log') # ax1 = plt.subplot(224) # plt.plot(data[:, 5], data[:, 8], 'o', markeredgecolor='g', markerfacecolor='none', label='EXP-PRE_M&H') plt.plot(x1, y1, '--', color='k') # , label='+30%') plt.plot(x1, y2, '--', color='k') # , label='-30%') # plt.legend(loc='best') # plt.xscale('log') # plt.yscale('log') # ax.set_aspect(1) # curr_time = datetime.datetime.now() # time_str=datetime.datetime.strftime(curr_time, '%Y%m%d%H%M') plt.savefig(time_str + 'res.pdf') plt.show() main()