一、问题背景结合最近一段时间学校布置的实验要求和数据本篇内容主要讲述如何通过Python用遗传算法对霍顿下渗方程的三个参数fc稳定下渗率f0初始下渗率k与土壤特性有关的经验常数并给出率定过程的动态图像展示。霍顿下渗方程其中f是t时刻的下渗率fc是稳定下渗率f0是初始下渗率t是时间k是常数。对上式进行积分可得下渗量的霍顿公式:实验思路首先写出目标函数及必要的约束条件运行遗传算法程序即可得到三个待定参数值。该目标函数是追求实际下渗量与模拟计算下渗量的差值的平方求和最小为最优目标目标函数如下式中Fi为实测下渗量fc、f0、k为霍顿下渗方程的三个参数。实验数据1t0已知代表土壤下渗率从初始下渗率f0对应的时刻到t时刻之间的时间长度t010min2需要说明的是,由于目标函数中的时间是t0t所以Fi应该直接选用W而不是W0。所以表1中的tmin、Wmm数据在后续需要用到。表1实验资料日期PmmRmmW0mmtminWmm序号12345611957.*.*23.42.710.0330.721959.6.123.42.59.5330.431959.7.486.617.810.92179.741959.8.514.41.412.9325.951959.8.745.45.89.21248.861960.9.323.61.211.0433.471962.7.268.52.68.11574二、代码展示import numpy as np import matplotlib.pyplot as plt from matplotlib import cm class GeneticAlgorithm: def __init__(self): self.DNA_SIZE 24 #DNA长度 self.POP_SIZE 80 #种群大小 self.CROSSOVER_RATE 0.6 #交叉率 self.MUTATION_RATE 0.01 #变异率 self.N_GENERATIONS 100 #迭代次数 self.X_BOUND [2, 3]# x的取值范围 self.Y_BOUND [0, 0.9]# y的取值范围 self.Z_BOUND [0.01, 0.3]# z的取值范围 self.sca None # 初始化 sca 属性 def F(self, x, y, z): f0 x fc y k z t0 10 E (30.7 - (fc * (t0 3) (1 / k) * (f0 - fc) * (1 - np.exp(-k * (t0 3))))) ** 2 \ (30.4 - (fc * (t0 3) (1 / k) * (f0 - fc) * (1 - np.exp(-k * (t0 3))))) ** 2 \ (79.7 - (fc * (t0 21) (1 / k) * (f0 - fc) * (1 - np.exp(-k * (t0 21))))) ** 2 \ (25.9 - (fc * (t0 3) (1 / k) * (f0 - fc) * (1 - np.exp(-k * (t0 3))))) ** 2 \ (48.8 - (fc * (t0 12) (1 / k) * (f0 - fc) * (1 - np.exp(-k * (t0 12))))) ** 2 \ (33.4 - (fc * (t0 4) (1 / k) * (f0 - fc) * (1 - np.exp(-k * (t0 4))))) ** 2 \ (74 - (fc * (t0 15) (1 / k) * (f0 - fc) * (1 - np.exp(-k * (t0 15))))) ** 2 # E (37.5 - (fc * (t0 5) (1 / k) * (f0 - fc) * (1 - np.exp(-k * (t0 5))))) ** 2 \ # (61.5 - (fc * (t0 12) (1 / k) * (f0 - fc) * (1 - np.exp(-k * (t0 12))))) ** 2 \ # (54.5 - (fc * (t0 8) (1 / k) * (f0 - fc) * (1 - np.exp(-k * (t0 8))))) ** 2 \ # (53.5 - (fc * (t0 9) (1 / k) * (f0 - fc) * (1 - np.exp(-k * (t0 9))))) ** 2 return E def get_fitness(self, pop): x, y, z self.translateDNA(pop) pred self.F(x, y, z) return np.max(pred) - pred 1e-3 def translateDNA(self, pop): x_pop pop[:, 0:self.DNA_SIZE] y_pop pop[:, self.DNA_SIZE:self.DNA_SIZE * 2] z_pop pop[:, self.DNA_SIZE * 2:] x x_pop.dot(2 ** np.arange(self.DNA_SIZE)[::-1]) / float(2 ** self.DNA_SIZE - 1) * (self.X_BOUND[1] - self.X_BOUND[0]) self.X_BOUND[0] y y_pop.dot(2 ** np.arange(self.DNA_SIZE)[::-1]) / float(2 ** self.DNA_SIZE - 1) * (self.Y_BOUND[1] - self.Y_BOUND[0]) self.Y_BOUND[0] z z_pop.dot(2 ** np.arange(self.DNA_SIZE)[::-1]) / float(2 ** self.DNA_SIZE - 1) * (self.Z_BOUND[1] - self.Z_BOUND[0]) self.Z_BOUND[0] return x, y, z def crossover_and_mutation(self, pop, CROSSOVER_RATE0.8): new_pop [] for father in pop: child father.copy() if np.random.rand() CROSSOVER_RATE: mother pop[np.random.randint(self.POP_SIZE)] cross_points np.random.randint(low0, highself.DNA_SIZE * 3) child[cross_points:] mother[cross_points:] self.mutation(child) new_pop.append(child) return new_pop def mutation(self, child, MUTATION_RATE0.01): for i in range(len(child)): if np.random.rand() MUTATION_RATE: mutate_point np.random.randint(0, self.DNA_SIZE * 3) child[mutate_point] child[mutate_point] ^ 1 def select(self, pop, fitness): idx np.random.choice(np.arange(self.POP_SIZE), sizeself.POP_SIZE, replaceTrue, pfitness / fitness.sum()) return pop[idx] def print_info(self, pop): fitness self.get_fitness(pop) max_fitness_index np.argmax(fitness) print(max_fitness:, fitness[max_fitness_index]) x, y, z self.translateDNA(pop) print(最优的基因型, pop[max_fitness_index]) print((x, y, z):, (x[max_fitness_index], y[max_fitness_index], z[max_fitness_index])) print(self.F(x[max_fitness_index], y[max_fitness_index], z[max_fitness_index])) def run(self): fig plt.figure() ax fig.add_subplot(111, projection3d) plt.ion() pop np.random.randint(2, size(self.POP_SIZE, self.DNA_SIZE * 3)) for _ in range(self.N_GENERATIONS): x, y, z self.translateDNA(pop) # Apply constraint -x y 0 valid_indices np.where(-x y 0) pop pop[valid_indices] x, y, z x[valid_indices], y[valid_indices], z[valid_indices] if len(pop) 0: print(No valid population left.) break if self.sca is not None:#实现动态展示 self.sca.remove() self.sca ax.scatter(x, y, z, cself.F(x, y, z), cmapcm.coolwarm) plt.show() plt.pause(0.1) pop np.array(self.crossover_and_mutation(pop, self.CROSSOVER_RATE)) fitness self.get_fitness(pop) pop self.select(pop, fitness) self.print_info(pop) plt.ioff() plt.show() def main(): ga GeneticAlgorithm() ga.run() if __name__ __main__: main()三、结果展示1迭代过程动态展示2参数率定结果