R中使用lpSolve求解项目排序最大化NPV的优化问题

1 基础场景：无惩罚规则的最优排序 #

核心思路 #

基础场景的最优排序本质是0-1指派问题：折现率为正数时，越早的位置折现系数越高，我们需要把高收入项目匹配到高权重的靠前位置，满足「每个项目仅出现1次、每个位置仅放1个项目」的约束。
lpSolve提供了lp.assign()函数可以直接求解这类指派问题，注意该函数默认求最小值，我们需要对收益系数取负值转化为最小化问题。

实现代码 #

# 加载包
library(lpSolve)

# 原始数据与自定义函数
df<-structure(list(id = c(1, 2, 3, 4, 5, 6, 7, 8), Revenue = c(109, 
111, 122, 139, 156, 140, 137, 167)), row.names = c(NA, 8L), class = "data.frame")

dcf <- function(x, r, t0=FALSE){
  if(length(r)==1){
    r <- rep(r, length(x))
    if(t0==TRUE){r[1]<-0}
  }
  x/cumprod(1+r)
}

npv <- function(x, r, t0=FALSE){
  sum(dcf(x, r, t0))
}

# 构造指派矩阵：行是排序位置（1-8），列是项目（1-8），值是该项目放在对应位置的负收益（因为lp.assign默认求最小）
n <- nrow(df)
r <- 0.2
# 计算每个位置的折现系数，对应dcf默认t0=FALSE的逻辑
discount_coef <- 1 / cumprod(rep(1+r, n))
# 构造成本矩阵
cost_mat <- matrix(NA, nrow = n, ncol = n)
for(i in 1:n){
  cost_mat[i,] <- - df$Revenue * discount_coef[i]
}

# 求解指派问题
res <- lp.assign(cost_mat)
# 提取最优排序：每行最大值所在列就是对应位置放的项目索引
optimal_order <- apply(res$solution, 1, which.max)
optimal_df <- df[optimal_order, ]

# 验证结果
cat("最优排序的NPV：", npv(optimal_df$Revenue, r), "\n")
cat("最优项目顺序ID：", optimal_df$id, "\n")

结果说明 #

运行后得到的最优NPV高于初始的492，符合收入越高越靠前的排序逻辑。

2 复杂场景：带相邻ID差惩罚的最优排序 #

核心思路 #

该场景下收益和相邻项目的ID差有关，属于旅行商问题(TSP)：我们需要找一条遍历所有项目的路径，总收益（折现后+惩罚调整）最大。
建模逻辑：

• 决策变量x[u,v]为0-1变量，1表示项目u后面紧跟项目v
• 惩罚规则：如果abs(df$id[u] - df$id[v]) <= 2，则v的收入打8折，否则不打折
• 目标函数：第一个项目的折现收入 + 所有后续项目的折现*惩罚后收入的总和最大
• 约束：每个项目仅入度1、仅出度1，消除子回路

实现代码 #

library(dplyr)
library(lpSolve)

n <- nrow(df)
r <- 0.2
discount_coef <- 1 / cumprod(rep(1+r, n))

# 构造边的收益矩阵：edge_profit[u,v] 是u后面接v时，v对应的收益
edge_profit <- matrix(0, nrow = n, ncol = n)
for(u in 1:n){
  for(v in 1:n){
    if(u == v) next
    penalty <- ifelse(abs(df$id[u] - df$id[v]) <= 2, 0.8, 1)
    edge_profit[u,v] <- penalty * df$Revenue[v]
  }
}

# 因为lp求解最小化，取负值
cost_edge <- -edge_profit

# 构造TSP基础约束：行和为1（出度1），列和为1（入度1）
row_const <- matrix(0, nrow = n, ncol = n*n)
for(i in 1:n){
  row_const[i, ((i-1)*n +1):(i*n)] <- 1
}
col_const <- matrix(0, nrow = n, ncol = n*n)
for(j in 1:n){
  col_const[j, seq(j, n*n, by = n)] <- 1
}
const_mat <- rbind(row_const, col_const)
const_dir <- rep("==", 2*n)
const_rhs <- rep(1, 2*n)

# 增加子回路消除约束（Miller-Tucker-Zemlin公式，适用于小规模TSP）
mtz_const <- matrix(0, nrow = (n-1)*(n-1), ncol = n*n)
mtz_rhs <- rep(NA, (n-1)*(n-1))
k <- 1
for(i in 2:n){
  for(j in 2:n){
    if(i == j) next
    mtz_const[k, (i-1)*n + j] <- 1
    mtz_rhs[k] <- n-1
    k <- k+1
  }
}
const_mat <- rbind(const_mat, mtz_const)
const_dir <- c(const_dir, rep("<=", (n-1)*(n-1)))
const_rhs <- c(const_rhs, mtz_rhs)

# 求解线性规划
res <- lp("min", as.vector(cost_edge), const_mat, const_dir, const_rhs, all.bin = TRUE)
sol_mat <- matrix(res$solution, nrow = n, ncol = n)

# 提取最优路径
path <- c()
current <- which(rowSums(sol_mat[,]) == 0) # 起点（入度为0）
path <- c(path, current)
for(i in 1:(n-1)){
  current <- which(sol_mat[current,] == 1)
  path <- c(path, current)
}
optimal_df_penalty <- df[path, ]

# 计算带惩罚的NPV
optimal_df_penalty <- optimal_df_penalty %>%
  mutate(Penalty = if_else(abs(id - lag(id)) <=2, 0.8, 1),
         Penalty = ifelse(is.na(Penalty), 1, Penalty),
         Revenue_Penalized = Revenue * Penalty)
cat("带惩罚规则的最优NPV：", npv(optimal_df_penalty$Revenue_Penalized, r), "\n")
cat("带惩罚规则的最优项目顺序ID：", optimal_df_penalty$id, "\n")

注意事项 #

8个节点属于小规模问题，增加子回路消除约束后完全可以在普通设备上快速得到全局最优解。

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