Python| 水文 |利用遗传算法对霍顿下渗方程进行参数率定

最新推荐文章于 2024-10-14 09:08:48 发布

摆烂老大

最新推荐文章于 2024-10-14 09:08:48 发布

阅读量925

点赞数 23

文章标签： python 算法数据分析

本文链接：https://blog.csdn.net/weixin_67437706/article/details/140421785

版权

一、问题背景

结合最近一段时间学校布置的实验要求和数据，本篇内容主要讲述如何通过Python用遗传算法对霍顿下渗方程的三个参数：fc（稳定下渗率），f0（初始下渗率），k（与土壤特性有关的经验常数），并给出率定过程的动态图像展示。

霍顿下渗方程：

其中，f是t时刻的下渗率，fc是稳定下渗率，f0是初始下渗率，t是时间，k是常数。对上式进行积分，可得下渗量的霍顿公式:

实验思路：首先写出目标函数，及必要的约束条件，运行遗传算法程序即可得到三个待定参数值。该目标函数是追求实际下渗量与模拟计算下渗量的差值的平方求和最小为最优目标

目标函数如下：

式中：Fi为实测下渗量；fc、f0、k为霍顿下渗方程的三个参数。

实验数据：

（1）t0已知，代表土壤下渗率从初始下渗率f0对应的时刻到t时刻之间的时间长度，t0=10min；

（2）需要说明的是,由于目标函数中的时间是t0+t，所以Fi应该直接选用W而不是W0。所以，表1中的t（min）、W（mm）数据在后续需要用到。

表1：实验资料

	日期	P（mm）	R（mm）	W0（mm）	t（min）	W（mm）
序号	（1）	（2）	（3）	（4）	（5）	（6）
1	1957..	23.4	2.7	10.0	3	30.7
2	1959.6.1	23.4	2.5	9.5	3	30.4
3	1959.7.4	86.6	17.8	10.9	21	79.7
4	1959.8.5	14.4	1.4	12.9	3	25.9
5	1959.8.7	45.4	5.8	9.2	12	48.8
6	1960.9.3	23.6	1.2	11.0	4	33.4
7	1962.7.2	68.5	2.6	8.1	15	74

二、代码展示

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_RATE=0.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(low=0, high=self.DNA_SIZE * 3)
                child[cross_points:] = mother[cross_points:]
            self.mutation(child)
            new_pop.append(child)
        return new_pop

    def mutation(self, child, MUTATION_RATE=0.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), size=self.POP_SIZE, replace=True, p=fitness / 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, projection='3d')
        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, c=self.F(x, y, z), cmap=cm.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()