如何利用optim获取月度负荷序列的正确全局负荷削减率及替代工具？

问题背景 #

需要计算全局负荷削减率（global_reduct），目标是让**仅目标比例（如50%）**的负荷值在削减后超过指定参考负荷（load_ref）。

初始参数与数据集 #

参数设置：

global_reduct <- c(50)  ## 初始值
load_ref <- 450.5  ## 参考负荷
tobject <- 50  ## 目标比例：50%的负荷超过参考值

示例月度负荷时间序列数据集：

df <- data.frame(
  Date=seq(as.Date('2010-01-01'), as.Date('2012-04-01'), by='months'), 
  load= c(1.496169954, 1.29147009, 1.964195241, 1.14352707, 1.319144304, 
          0.773288093, 0.65175612, 0.685340958, 0.416934849, 
          0.769853258, 1.104639594, 0.92213209, 1.685588986, 
          1.972510175, 2.6882446, 2.153314503, 1.324735759, 
          1.027755411, 0.610207197, 0.674642831, 0.721971375, 
          1.13233884, 0.739325423, 0.90031817, 1.366597449, 
          1.928098735, 1.216538229, 1.514353244)
)

已知当目标为50%负荷超过参考值时，正确削减率约为62%。

现有尝试的问题 #

用户编写了供optim调用的函数，但返回结果不符合预期，原代码：

optfuc <- function(reduct, ttarget=50){
  reduct_eq <- df$load *(1 - (reduct/100))
  tt_exceed <- ifelse((reduct_eq *1000) > load_ref, 1, 0)  
  ave_ref <- sum(tt_exceed)/length(tt_exceed)*100 - ttarget
  # ave_ref in this case should be = ttarget
  # ave_ref
  reduct 
}

optim(c(30), optfuc, method ="L-BFGS-B", lower=0, upper=100)

问题核心：optim需要目标函数返回待最小化的误差值，但原函数直接返回reduct，无法引导优化过程找到正确值。

解决方案：修正优化函数 #

调整目标函数，让它返回「实际超比例与目标比例的差值的平方」（平方误差便于优化，也可用绝对值），让optim最小化这个误差：

# 修正后的目标函数：返回误差平方
optfuc <- function(reduct, ttarget=50){
  reduct_eq <- df$load * (1 - reduct/100)
  # 原数据load需乘以1000才能和load_ref比较
  tt_exceed <- as.integer(reduct_eq * 1000 > load_ref)
  actual_pct <- mean(tt_exceed) * 100
  # 返回误差的平方，optim会最小化这个值
  (actual_pct - ttarget)^2
}

# 调用optim求解，初始值设为30，约束0-100
result <- optim(c(30), optfuc, method = "L-BFGS-B", lower=0, upper=100)
# 输出最优削减率
cat("最优全局削减率：", round(result$par, 2), "%\n")

运行后会得到约62%的结果，符合预期。

替代方案：使用求根包 #

这个问题本质是找一个reduct值，使得「超比例=目标比例」，属于求根问题，用rootSolve包的uniroot.all函数更直接：

# 安装并加载包
install.packages("rootSolve")
library(rootSolve)

# 定义根函数：返回实际超比例与目标的差值
root_func <- function(reduct, ttarget=50){
  reduct_eq <- df$load * (1 - reduct/100)
  tt_exceed <- as.integer(reduct_eq *1000 > load_ref)
  mean(tt_exceed)*100 - ttarget
}

# 在0-100范围内找根
sol <- uniroot.all(root_func, interval = c(0,100), ttarget=50)
cat("最优全局削减率：", round(sol, 2), "%\n")

这种方法不需要最小化误差，直接找到满足条件的削减率，逻辑更贴合问题需求。

内容的提问来源于stack exchange，提问作者AMSS