# 清理环境（可选）
rm(list = ls())
# 加载必要的包
library(pso)
### --- 需要手动输入的部分开始 --- ###
# 1. 设置工作目录（包含DSSAT运行文件的位置）
setwd('C:/DSSAT48/Maize/') # ⚠️ 修改为您的实际路径
# 2. 定义文件路径（根据您的DSSAT安装位置修改）
DSSATcul_Path <- 'C:/DSSAT48/Genotype/MZCER048.CUL' # ️ 修改为实际CUL文件路径
DSSATeco_Path <- 'C:\\DSSAT48\\Genotype\\MZCER048.ECO' # ️ 修改为实际ECO文件路径
soil_file_Path <- 'C:\\DSSAT48\\Soil\\SOIL.SOL' # ️ 修改为实际土壤文件路径
bat_file_path <- "C:/DSSAT48/Maize/批处理.bat" # ️ 修改为实际批处理文件路径
# 3. 输入观测数据（根据您的实验数据修改）
M_DAS <- as.numeric(c(95,95,95,95,95,95,95,95,95))
HWAMM <- as.numeric(c(4125,5825,6571,5676,5336,6181,4042,3819,4616))
CWAMM <- as.numeric(c(13185,16951,19518,13269,16830,16686,11776,14022,13695))
# 4. 定义参数范围（6个参数，每行一个参数的上下界）
param_ranges <- matrix(c(
5.0, 450.0, # ️ 参数1范围
0.000, 2.000, # ️ 参数2范围
580.0, 999.0, # ️ 参数3范围
248.0, 990.0, # ️ 参数4范围
5.00, 16.50, # ⚠️ 参数5范围
38.00, 75.00 # ⚠️ 参数6范围
), ncol = 2, byrow = TRUE)
### --- 需要手动输入的部分结束 --- ###
### 以下为自动处理部分（一般不需要修改）###
# 预计算常量值（观测数据的均值）
mean_M_DAS <- mean(M_DAS)
mean_HWAMM <- mean(HWAMM)
mean_CWAMM <- mean(CWAMM)
# 读取CUL文件并提取特定基因型行
cul_file <- readLines(DSSATcul_Path)
Genename <- '999996'
culNo0 <- grep(pattern = Genename, cul_file)
cullines <- substr(cul_file[culNo0], 1, 36) # 保留前36个字符
# 预先读取ECO和土壤文件（仅读取，后续不再修改）
eco_file <- readLines(DSSATeco_Path)
sol_file <- readLines(soil_file_Path)
# 设置参数边界
lower_bounds <- param_ranges[, 1]
upper_bounds <- param_ranges[, 2]
# 初始化计数器
iteration_counter <- 0
# 定义运行DSSAT模型的函数（优化版）
rundssat_opt <- function(params) {
# 参数边界处理
params <- pmin(pmax(params, lower_bounds), upper_bounds)
# 精确复制原始代码的参数格式化逻辑
formatted_cul_params <- sapply(params, function(x) {
param_str <- format(x, scientific = FALSE, nsmall = 4) # 确保四位小数
substr(param_str, 1, 5) # 取前五位
})
# 拼接参数字符串
new_cul_params_str <- paste(formatted_cul_params, collapse = " ")
# 更新CUL文件内容（精确复制原始代码）
new_cul_lines <- paste(cullines, new_cul_params_str, sep = " ")
cul_file[culNo0] <- new_cul_lines
# 写入CUL文件
writeLines(cul_file, con = DSSATcul_Path)
# 运行DSSAT模型（批处理文件）
system(bat_file_path, intern = TRUE, show.output.on.console = FALSE)
# 检查结果文件是否存在
eval_path <- "C:/DSSAT48/Maize/Evaluate.OUT"
if (!file.exists(eval_path)) {
return(list(M_DATE = NA, HWAMS = NA, CWAMS = NA))
}
# 读取结果文件（只读前n行，然后取第 [到] 行）
eval_data <- suppressWarnings(readLines(eval_path, n = 12)[4:12])
# 辅助函数：提取指定位置的数值
extract_cols <- function(data, start, end) {
as.numeric(trimws(substr(data, start, end)))
}
# 提取模拟值
M_DATE <- extract_cols(eval_data, 81, 83)
HWAMS <- extract_cols(eval_data, 95, 99)
CWAMS <- extract_cols(eval_data, 175, 179)
return(list(M_DATE = M_DATE, HWAMS = HWAMS, CWAMS = CWAMS))
}
# 定义适应度函数（优化版）
fitness_function_opt <- function(params) {
iteration_counter <<- iteration_counter + 1
cat("Iteration:", iteration_counter, "\n")
# 运行DSSAT模型获取模拟值
sim_values <- rundssat_opt(params)
# 如果模拟失败，返回一个大值（表示差）
if (any(is.na(unlist(sim_values)))) {
cat("Simulation failed. Assigning high penalty.\n")
return(1e6)
}
# 计算均方误差
sq_diff_M_DATE <- (sim_values$M_DATE - M_DAS)^2
sq_diff_HWAMS <- (sim_values$HWAMS - HWAMM)^2
sq_diff_CWAMS <- (sim_values$CWAMS - CWAMM)^2
# 计算RMSE
rmse_m_date <- sqrt(mean(sq_diff_M_DATE))
rmse_hwams <- sqrt(mean(sq_diff_HWAMS))
rmse_cwams <- sqrt(mean(sq_diff_CWAMS))
# 计算NRMSE（归一化均方根误差）
nrmse_m_date <- rmse_m_date / mean_M_DAS
nrmse_hwams <- rmse_hwams / mean_HWAMM
nrmse_cwams <- rmse_cwams / mean_CWAMM
# 加权NRMSE
weights <- c(0.2, 0.4, 0.3) # 权重分配
weighted_nrmse <- sum(weights * c(nrmse_m_date, nrmse_hwams, nrmse_cwams))
# 计算R²（决定系数）仅针对HWAMS
if (sd(HWAMM) > 0 && sd(sim_values$HWAMS) > 0) {
r_squared <- cor(sim_values$HWAMS, HWAMM)^2
# 适应度函数：加权NRMSE减去R²
combined_fitness <- 0.6 * weighted_nrmse - 0.4 * r_squared
} else {
# 如果无法计算R²，则仅使用加权NRMSE
combined_fitness <- 0.6 * weighted_nrmse + 0.4 * 1e3 # 惩罚项
}
cat("Fitness:", combined_fitness,
"| Weighted NRMSE:", weighted_nrmse,
"| R²:", ifelse(exists("r_squared"), r_squared, "NA"), "\n")
return(combined_fitness)
}
### PSO优化部分 ###
# 设置PSO参数（可调整）
pso_params <- list(
maxit = 200, # 最大迭代次数（⚠️ 可根据需要调整）
#s = 30, # 种群大小（粒子数）
#reltol = 1e-6, #相对收敛
#abstol = 1e-10, # 收敛阈值（当适应度变化小于此值时停止）
trace = 1, # 显示迭代过程
REPORT = 5 # 每5次迭代报告一次
)
# 重置计数器
iteration_counter <- 0
# 初始化粒子群（使用参数范围的中点作为初始值）
initial_par <- apply(param_ranges, 1, mean) # 取每个参数范围的中点
# 运行PSO优化
pso_result <- psoptim(
par = initial_par, # 初始参数值
fn = fitness_function_opt,
lower = lower_bounds,
upper = upper_bounds,
control = pso_params
)
### 结果输出部分 ###
# 输出优化结果
cat("\n===== Optimization Results =====\n")
cat("Optimization finished after", iteration_counter, "iterations\n")
cat("Best fitness value:", pso_result$value, "\n\n")
# 格式化输出最佳参数（精确复制原始格式化方式）
best_params_formatted <- sapply(pso_result$par, function(x) {
param_str <- format(x, scientific = FALSE, nsmall = 4)
substr(param_str, 1, 5)
})
cat("Optimized parameters (formatted):\n")
print(best_params_formatted)
# 使用最佳参数运行一次模型
cat("\nRunning simulation with optimized parameters...\n")
final_sim <- rundssat_opt(pso_result$par)
# 输出模拟结果与观测值的对比
cat("\n===== Simulation vs Observation =====\n")
cat("M_DATE (Simulated):", final_sim$M_DATE, "\n")
cat("M_DAS (Observed):", M_DAS, "\n\n")
cat("HWAMS (Simulated):", final_sim$HWAMS, "\n")
cat("HWAMM (Observed):", HWAMM, "\n\n")
cat("CWAMS (Simulated):", final_sim$CWAMS, "\n")
cat("CWAMM (Observed):", CWAMM, "\n")
你学习下这个代码，然后这个代码命名为“代码一”
本文介绍了C语言中的常量种类，包括字面常量、const修饰的常变量和枚举常量，以及注释的两种方式。同时，讲解了转义字符的用法，如
用于换行，并通过printf函数示例展示了如何打印不同类型的值，如%d用于整数，%c用于字符，%s用于字符串。还提到了八进制表示字符的方法及其与ASCII码的关系。
扫一扫
532
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
