# -*- coding: utf-8 -*-
import numpy as np
import os
import pandas as pd
from deap import base, creator, tools, algorithms
import subprocess
import random
import openpyxl
from netCDF4 import Dataset
import multiprocessing
import configparser
import shutil
global result_list
result_list = {}
excel_template = r'E:\cwat_new\diversion_optimization\diversion_template.xlsx'
ini_file = r'E:\cwat_new\diversion_optimization\optim.ini'
total_water = 15128
used_ids = set()
def generate_unique_id():
unique_id = random.randint(1, 100000000)
while unique_id in used_ids:
unique_id = random.randint(1, 100000000)
used_ids.add(unique_id)
return unique_id
def read_initial_values(excel_file_path, column_name):
workbook = openpyxl.load_workbook(excel_file_path)
sheet = workbook.active
# 查找指定列的索引
column_index = None
for cell in sheet[1]:
if cell.value == column_name:
column_index = cell.column_letter
break
if column_index is None:
print(f"找不到名为 '{column_name}' 的列标题")
return []
# 获取指定列下的所有值
column_values = []
for cell in sheet[column_index][1:]:
column_values.append(cell.value)
return column_values
def evaluate_model_output(nc_file):
with Dataset(nc_file, 'r') as nc:
times = nc.variables['time'][:]
data_var = nc.variables['unmetDemandM3_annualavg']
# 遍历每个时间段，计算总和
for i, _ in enumerate(times):
data_at_time = data_var[i, :, :]
data_at_time = data_at_time * 365
# 计算当前时间段的总和
total = np.sum(data_at_time)
#print(f"Time {i} total: {total}")
total = tuple([total])
return total
# 更新Excel文件中的月份分配数据
def update_excel_data(filename, allocations):
df = pd.DataFrame(allocations, columns=['Allocation'])
df['Month'] = range(1, 13)
with pd.ExcelWriter(filename, engine='openpyxl', mode='w') as writer:
df.to_excel(writer, index=False)
def mutate_individual(individual, mutation_probability):
for i in range(len(individual)):
if random.random() < mutation_probability:
individual[i] += random.gauss(0, 100) # 变异操作
# 修正子代，确保在非负范围内
individual[i] = max(individual[i], 0)
return [individual]
def cxBlendBounded(ind1, ind2, alpha=0.5, low=0):
"""执行cxBlend交叉操作，并确保子代的值位于[low, up]范围内。"""
for i in range(len(ind1)):
gamma = (1. + 2. * alpha) * random.random() - alpha
ind1[i] = (1. - gamma) * ind1[i] + gamma * ind2[i]
ind2[i] = gamma * ind1[i] + (1. - gamma) * ind2[i]
# 确保值不小于下限
ind1[i] = max(low, ind1[i])
ind2[i] = max(low, ind2[i])
return ind1, ind2
# 运行水文模型的命令行函数
def run_hydrological_model(unique_id, monthly_allocation):
folder_name = f'E:\\cwat_new\\diversion_optimization\\{unique_id}'
os.makedirs(folder_name, exist_ok=True)
output_folder = os.path.join(folder_name, 'output')
os.makedirs(output_folder, exist_ok=True)
excel_file = os.path.join(folder_name, f'diversion_{unique_id}.xlsx')
# 使用模板文件创建新的Excel
shutil.copy(excel_template, excel_file)
#创建新的ini文件
config = configparser.ConfigParser()
config.optionxform = lambda option: option
config.read(ini_file)
# 修改参数值
config.set('WATERDIVERSION', 'monthly_water_quota', excel_file)
config.set('FILE_PATHS', 'PathOut', output_folder)
# 保存修改后的ini文件
ini_fileC = os.path.join(folder_name, f'optim_{unique_id}.ini')
with open(ini_fileC, 'w') as configfile:
config.write(configfile)
# 更新Excel文件以供模型使用
update_excel_data(excel_file, monthly_allocation)
# 运行水文模型，假设模型的可执行文件名为 "hydro_model"，并且它使用excel文件
subprocess.run(['python', r'E:\cwat_new\CWatM-main\run_cwatm.py', ini_fileC], check=True)
nc_file = os.path.join(output_folder, 'unmetDemandM3_annualavg.nc')
result_list[int(evaluate_model_output(nc_file)[0])] = unique_id
# 从生成的nc文件中读取评价指标
return evaluate_model_output(nc_file), unique_id
# 评价函数，目标是最小化评价指标
def evaluate(individual):
unique_id = generate_unique_id()
# 确保分配水量不超过限制
if not isinstance(individual, list):
individual = list(individual)
# 确保 individual 中的元素都是整数
individual = [int(x) for x in individual]
# 计算分配方案的总和
allocation_sum = sum(individual)
penalty = abs(allocation_sum - total_water) * 1000000
if allocation_sum > total_water:
return 1e10, # 返回一个很大的值，表示不可行的解
result, unique_id = run_hydrological_model(unique_id, individual)
return result[0] + penalty,
# 设置遗传算法
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individual", list, fitness=creator.FitnessMin)
toolbox = base.Toolbox()
toolbox.register("attr_float", random.uniform, 0, total_water/12) # 假设平均每月分配
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_float, n=12)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
toolbox.register("evaluate", evaluate)
toolbox.register("mate", tools.cxBlend, alpha=0.5)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=100, indpb=0.2)
toolbox.register("select", tools.selTournament, tournsize=3)
# 遗传算法参数
population_size = 1
crossover_probability = 0.7
mutation_probability = 0.2
number_of_generations = 0
toolbox.register("mutate", mutate_individual, mutation_probability=mutation_probability)
toolbox.register("mate", cxBlendBounded, alpha=0.5, low=0)
if __name__ == "__main__":
manager = multiprocessing.Manager()
result_list = manager.dict()
pool_size = int(multiprocessing.cpu_count() * 0.8)
pool = multiprocessing.Pool(processes=pool_size)
toolbox.register("map", pool.map)
population = toolbox.population(n=population_size)
#设置初始值
initial_values = read_initial_values(excel_template, 'Allocation')
population[0][:] = initial_values
final_population, logbook = algorithms.eaSimple(population, toolbox, cxpb=crossover_probability, mutpb=mutation_probability, ngen=number_of_generations, verbose=True)
# 找到最优解
best_ind = tools.selBest(population, 1)[0]
best_fitness = best_ind.fitness.values[0]
best_unique_id = result_list[int(best_fitness)]
print("Best Individual is: ", best_ind)
print("Best Individual fitness:", best_fitness)
print("Best Individual ID is:", best_unique_id)
30 鱼币
( 粤ICP备18085999号-1 | 粤公网安备 44051102000585号)
狗仔卡
发表于 2024-4-2 07:59:53
显身卡
发表于 2024-4-2 18:35:52
