如何用R实现满足行列和约束的n×n矩阵线性规划优化碳排放

• 需求：在R语言中开发程序，优化243个国家间航空、水运、公路、铁路4种运输方式的贸易布局，以最小化碳排放。具体要求为选择4个贸易矩阵，使矩阵元素与对应排放系数逐元素相乘的总和最小，且总贸易的行列和需与样本一致。
• 现有问题：编写的3×3规模示例代码无法运行，尝试通过求和矩阵设置行列和约束未成功，最终需将方案扩展至243×243规模，请求排查问题并提供解决方案。

原示例代码 #

# Test Emission by MOT 1 Matrix
E1 <- matrix(c(0,1,1,1,0,1,1,1,0), nrow = 3)

print(E1)

# Test Emission by MOT 2 Matrix
E2 <- matrix(c(0,2,2,2,0,2,2,2,0), nrow = 3)

print(E2)

# Test Distance between three countries matrix
D <- matrix(c(0,1,2,1,0,1,2,1,0), nrow = 3)

print(D)


objective.in <- function(values, A, B) {
  X <- matrix(values[1:9], nrow = 3, ncol = 3)
  Y <- matrix(values[10:18], nrow = 3, ncol = 3)
  return(sum(E1 * D * X + E2 * D * Y))
}

Trade_MOT1_Hypothetical <- matrix(c(0,4,6,4,0,5,2,1,0), nrow = 3)
Trade_MOT2_Hypothetical <- matrix(c(0,8,9,1,0,17,3,1,0), nrow = 3)

print(Trade_MOT1_Hypothetical)
print(Trade_MOT2_Hypothetical)

test_input <- c(Trade_MOT1_Hypothetical, Trade_MOT2_Hypothetical)

print(test_input)

print(objective.in(test_input))

# Summation matrix for the first column sum constraint

c1 <- matrix(c(1, 1, 1, 0, 0, 0, 0, 0, 0), nrow = 9)

# Summation matrix for the second column sum constraint
c2 <- matrix(c(0, 0, 0, 1, 1, 1, 0, 0, 0), nrow = 9)

# Summation matrix for the third column sum constraint
c3 <- matrix(c(0, 0, 0, 0, 0, 0, 1, 1, 1), nrow = 9)

# Summation matrix for the first row sum constraint
r1 <- matrix(c(1, 0, 0, 1, 0, 0, 1, 0, 0), nrow = 9)

# Summation matrix for the second row sum constraint
r2 <- matrix(c(0, 1, 0, 0, 1, 0, 0, 1, 0), nrow = 9)

# Summation matrix for the third row sum constraint
r3 <- matrix(c(0, 0, 1, 0, 0, 1, 0, 0, 1), nrow = 9)

const.mat <- matrix(c(r1, r2, r3, c1, c2, c3, r1, r2, r3, c1, c2, c3), nrow = 6)

print(const.mat)

const.dir <- c("=", "=", "=", "=", "=", "=")

const.rhs <- c(1000, 1000, 1000, 2000, 2000, 2000)

optimum <-  lp(direction = "min", objective.in, const.mat, const.dir, const.rhs)

optimum$solution

print(optimum$solution)

1. 目标函数参数不匹配：lp函数要求目标函数仅接受决策变量向量作为唯一参数，但原代码中objective.in额外定义了A, B参数，导致调用时参数错误。
2. 约束矩阵构造错误：原代码将多个9×9矩阵拼接成6行矩阵，维度完全不符合要求。正确的约束矩阵应为每行对应一个约束，每列对应一个决策变量的结构，3×3两种运输方式的场景下，约束矩阵应为12行（6个总贸易行列和约束）×18列（18个决策变量）。
3. 约束逻辑偏离需求：原代码设置单个运输方式的行列和为固定值，但实际需求是所有运输方式的贸易矩阵之和的行列和与样本一致，约束逻辑完全错误。
4. 缺少依赖包调用：代码使用lp函数但未导入lpSolve包，会直接报错。

1. 可运行的3×3规模代码 #

# 导入依赖包
library(lpSolve)

# 排放系数矩阵（两种运输方式）
E1 <- matrix(c(0,1,1,1,0,1,1,1,0), nrow = 3)
E2 <- matrix(c(0,2,2,2,0,2,2,2,0), nrow = 3)
# 国家间距离矩阵
D <- matrix(c(0,1,2,1,0,1,2,1,0), nrow = 3)

# 目标函数：仅接受决策变量向量作为参数
objective.in <- function(values) {
  n <- 3
  # 提取两种运输方式的贸易矩阵
  X <- matrix(values[1:(n^2)], nrow = n, ncol = n)
  Y <- matrix(values[(n^2+1):(2*n^2)], nrow = n, ncol = n)
  # 计算总碳排放
  return(sum(E1 * D * X + E2 * D * Y))
}

# 样本总贸易的行列和（所有运输方式之和的行列约束）
sample_row_sums <- c(10, 9, 3)
sample_col_sums <- c(6, 5, 11)

n <- 3
num_vars <- 2 * n^2  # 决策变量总数：2个3×3矩阵共18个变量

# 构造约束矩阵：每行对应一个约束，每列对应一个决策变量
const_mat <- matrix(0, nrow = 2*n, ncol = num_vars)

# 行和约束：X+Y的行和等于样本行和
for (i in 1:n) {
  # 标记X矩阵第i行的变量位置
  const_mat[i, (i-1)*n + 1:i*n] <- 1
  # 标记Y矩阵第i行的变量位置
  const_mat[i, n^2 + (i-1)*n +1:i*n] <- 1
}

# 列和约束：X+Y的列和等于样本列和
for (j in 1:n) {
  row_idx <- n + j
  # 标记X矩阵第j列的变量位置
  const_mat[row_idx, seq(j, n^2, by = n)] <- 1
  # 标记Y矩阵第j列的变量位置
  const_mat[row_idx, n^2 + seq(j, n^2, by = n)] <- 1
}

# 约束方向和右侧值
const_dir <- rep("=", 2*n)
const_rhs <- c(sample_row_sums, sample_col_sums)

# 贸易量非负约束
lower_bounds <- rep(0, num_vars)

# 求解线性规划
optimum <- lp(direction = "min", 
              objective.in = objective.in, 
              const.mat = const_mat, 
              const.dir = const_dir, 
              const.rhs = const_rhs,
              lower.tail = lower_bounds)

# 输出结果
cat("优化状态：", optimum$status, "\n")
if (optimum$status == 0) {
  cat("最小碳排放：", optimum$objval, "\n")
  # 提取优化后的贸易矩阵
  X_opt <- matrix(optimum$solution[1:n^2], nrow = n, ncol = n)
  Y_opt <- matrix(optimum$solution[(n^2+1):(2*n^2)], nrow = n, ncol = n)
  cat("运输方式1的优化贸易矩阵：\n")
  print(X_opt)
  cat("运输方式2的优化贸易矩阵：\n")
  print(Y_opt)
}

2. 扩展到243×243规模的注意事项 #

• 求解器选择：243个国家×4种运输方式的场景下，决策变量数量达4*243^2 = 235224，lpSolve无法处理该规模，建议使用ROI包配合glpk/cplex后端，或lpsolveAPI调用底层API。
• 内存优化：使用Matrix包的稀疏矩阵格式存储约束矩阵，避免完整矩阵占用大量内存。
• 变量简化：设置对角线元素（国家内部贸易）的上下界为0，减少有效变量数量。
• 并行加速：启用求解器的并行计算功能，提升求解效率。

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