如何在R中提升可再生能源电站电池系统顺序分析效率

问题背景 #

为可再生能源电站设计电池系统，满足机器1000kWh左右的供电需求：

• 发电量超过需求时，剩余电量用于给100kW电池充电
• 发电量不足且电池有剩余容量时，由电池补能
• 现有数据集：1040条路径、8760小时的energia矩阵，需生成carga_bateria（电池电量）和更新后的energia数据集
• 当前实现：逐小时for循环+多条件判断，因顺序依赖无法并行，模拟效率偏低

现有代码的性能瓶颈 #

1. 使用data.frame存储大规模数值矩阵（1040×8760），R中data.frame的元素访问/修改效率远低于原生矩阵
2. 循环中重复计算carga_bateria[,i-1]*(1-0.02/100)等中间变量，造成不必要的运算浪费
3. 条件判断的逻辑分支存在冗余，逐行赋值的验证统计操作效率低

优化方案 #

1. 改用矩阵存储替代Data Frame #

R的原生矩阵是连续内存存储的同类型数据，访问和修改速度远快于data.frame，直接初始化矩阵：

carga_bateria <- matrix(0, nrow = 1040, ncol = 8760)
validacion <- matrix(0, nrow = 5, ncol = 8760) # 预分配验证矩阵内存

2. 预计算重复中间变量 #

将循环中多次用到的损耗后电量、可用放电量等变量提前计算，避免重复运算：

for (i in 2:8760) {  
  # 预计算前一小时电池扣除损耗后的电量
  prev_bat <- carga_bateria[,i-1] * (1 - 0.02/100)
  # 预计算电池可用放电量（扣除20%最低容量）
  available_discharge <- pmax(prev_bat - n_bat*7.55*0.2, 0)
  
  # 条件1：发电量超需求，剩余电量充电
  cond1 <- energia[,i] > pot_elec
  carga_bateria[cond1,i] <- prev_bat[cond1] + pmin(n_bat*7.55 - prev_bat[cond1], (energia[cond1,i] - pot_elec)*0.85)
  
  # 条件2：发电量+可用放电量仍不足需求30%，优先充电
  cond2 <- (energia[,i] + available_discharge) < pot_elec*0.3
  carga_bateria[cond2,i] <- prev_bat[cond2] + pmin(n_bat*7.55 - prev_bat[cond2], energia[cond2,i]*0.85)
  
  # 条件3：电量缺口可由电池完全覆盖，部分放电补能
  cond3 <- pot_elec > energia[,i] & (pot_elec - energia[,i]) < available_discharge
  carga_bateria[cond3,i] <- prev_bat[cond3] - (pot_elec - energia[cond3,i])
  energia[cond3,i] <- pot_elec
  
  # 条件4：发电量+可用放电量介于需求30%-100%，电池完全放电补能
  cond4 <- (energia[,i] + available_discharge) > pot_elec*0.3 & (energia[,i] + available_discharge) < pot_elec
  carga_bateria[cond4,i] <- n_bat*7.55*0.2
  energia[cond4,i] <- energia[cond4,i] + (prev_bat[cond4] - n_bat*7.55*0.2)
  
  # 其他场景：仅扣除电池损耗
  cond99 <- !cond1 & !cond2 & !cond3 & !cond4
  carga_bateria[cond99,i] <- prev_bat[cond99]
  
  # 批量赋值验证统计，替代逐行赋值
  validacion[,i] <- c(sum(cond1), sum(cond2), sum(cond3), sum(cond4), sum(cond1|cond2|cond3|cond4))
}

3. 用Rcpp实现核心循环（终极优化） #

如果上述优化仍无法满足速度需求，可将核心循环逻辑用C++实现，利用编译型语言的效率优势：

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
List simulate_battery(NumericMatrix energia, double pot_elec, double n_bat) {
  int n_paths = energia.nrow();
  int n_hours = energia.ncol();
  NumericMatrix carga_bateria(n_paths, n_hours);
  NumericMatrix new_energia = clone(energia);
  
  const double bat_cap = n_bat * 7.55;
  const double min_bat = bat_cap * 0.2;
  const double loss_factor = 0.9998; // 1 - 0.02/100
  const double charge_eff = 0.85;

  for (int i = 1; i < n_hours; ++i) {
    for (int j = 0; j < n_paths; ++j) {
      const double prev_bat = carga_bateria(j, i-1) * loss_factor;
      const double available_discharge = std::max(prev_bat - min_bat, 0.0);
      double curr_energia = new_energia(j, i);

      if (curr_energia > pot_elec) {
        const double charge_amount = std::min(bat_cap - prev_bat, (curr_energia - pot_elec)*charge_eff);
        carga_bateria(j, i) = prev_bat + charge_amount;
      } else if (curr_energia + available_discharge < pot_elec * 0.3) {
        const double charge_amount = std::min(bat_cap - prev_bat, curr_energia * charge_eff);
        carga_bateria(j, i) = prev_bat + charge_amount;
      } else if (pot_elec > curr_energia && (pot_elec - curr_energia) < available_discharge) {
        carga_bateria(j, i) = prev_bat - (pot_elec - curr_energia);
        new_energia(j, i) = pot_elec;
      } else if (curr_energia + available_discharge > pot_elec * 0.3 && curr_energia + available_discharge < pot_elec) {
        carga_bateria(j, i) = min_bat;
        new_energia(j, i) = curr_energia + (prev_bat - min_bat);
      } else {
        carga_bateria(j, i) = prev_bat;
      }
    }
  }

  return List::create(Named("carga_bateria") = carga_bateria, Named("energia") = new_energia);
}

在R中调用方式：

library(Rcpp)
sourceCpp("battery_sim.cpp")

# 运行模拟
result <- simulate_battery(energia, pot_elec, n_bat)
carga_bateria <- result$carga_bateria
energia <- result$energia

4. 内存优化细节 #

• 避免动态扩展对象，提前预分配carga_bateria、validacion等矩阵的内存
• 尽量减少数据复制操作，比如直接修改矩阵元素而非创建新对象

内容的提问来源于stack exchange，提问作者Noe Sebastián Medina Muñoz