为何Python混合整数线性规划求解速度远慢于R？求优化方案

我有一个需运行数千次的大规模*混合整数线性规划（MILP）*问题，速度是核心需求。该问题需同时在R和Python中实现，其中R版代码求解速度极快，仅需约0.002秒；但Python版求解耗时在0.05至0.8秒之间，有时甚至更久。

我绝不认为R的效率比Python高25倍以上，毕竟据我所知两者均基于C代码实现求解器。现寻求帮助优化Python代码以提升速度，推测需更换当前使用的包，目标是让Python代码稳定在0.01秒内完成求解。

以下是求解同一问题的R和Python示例代码（以构建符合约束的梦幻篮球阵容为场景：500名球员中选8人，每个位置至少1人、最多3人（部分球员为双位置），总薪资不超50000；随机生成位置/薪资数据，重点关注求解速度而非实际解）：

R 代码 #

library(lpSolve)
library(dplyr)

# Make dummy dataframe with player positions, salaries, and projected points
Positions <- c('PG', 'SG', 'SF', 'PF', 'C', 'PG/SG', 'SG/SF', 'PF/C')
Position <- sample(Positions, 500, replace = TRUE)
Salary <- sample(seq(3000, 12000, by=100), 500, replace=TRUE)
playerTable <- data.frame(ID = 1:500, Position, Salary) %>% 
  mutate(ProjPts = round(Salary/1000*5))

# Make positional constraint vectors
ConVec_Position <- playerTable %>% 
  select(Position) %>% 
  mutate(AnyPG =  grepl('PG', Position)*1,
         AnySG = grepl('SG', Position)*1,
         AnySF = grepl('SF', Position)*1,
         AnyPF = grepl('PF', Position)*1,
         AnyC = grepl('C', Position)*1,
         
         orPGSGSFPFC = (AnyPG + AnySG + AnySF + AnyPF + AnyC > 0)*1
  ) %>% 
  select(-Position)

# Make salary constraint vectors
ConVec_Salary <- playerTable$Salary

# Get Objective (projected points)
ProjectedPoints <- playerTable$ProjPts

# Make constraint values and directions
ConValGT <- c(1,1,1,1,1, 8) # Must have at least one of each position, and 8 total
ConValLT <- c(3,3,3,3,2, 8, 50000) # Can't have more than 3 of each position, 8 total, and 50000 salary
ConDirGT <- c(rep(">=", ncol(ConVec_Position))) # Make >= directions for each position constraint
ConDirLT <- c(rep("<=", ncol(ConVec_Position) + 1)) # Make <= directions for each position constraint, plus 1 for Salary

# Compile constraint vectors, values, and directions
ConVec_Final <- cbind(ConVec_Position, ConVec_Position, ConVec_Salary)
ConDir_Final <- c(ConDirGT, ConDirLT)
ConVal_Final <- c(ConValGT, ConValLT)

start_time <- proc.time()

sol <- lp("max",
          objective.in = ProjectedPoints,
          const.mat = t(ConVec_Final),
          const.dir = ConDir_Final,
          const.rhs = ConVal_Final,
          binary.vec = 1:length(ProjectedPoints) # All decisions are binary
)

cat("\n--- It took ", round((proc.time() - start_time)[3], 4), " seconds to optimize lineups ---
", sep = "")

Python 原代码 #

import pandas as pd
import numpy as np
from lp_solve import *
import time

# Make dummy dataframe with player positions, salaries, and projected points
ID = range(500)
Positions = ['PG', 'SG', 'SF', 'PF', 'C', 'PG/SG', 'SG/SF', 'PF/C']
Position = np.random.choice(Positions,500)
Salary = np.random.choice(range(3000,12000,100),500)
playerTable = pd.DataFrame({'ID': ID, 'Position': Position, 'Salary': Salary})
playerTable['ProjPts'] = playerTable['Salary']/1000*5

# Make positional constraint vectors
AnyPG = playerTable.Position.str.contains('PG')*1
AnySG = playerTable.Position.str.contains('SG')*1
AnySF = playerTable.Position.str.contains('SF')*1
AnyPF = playerTable.Position.str.contains('PF')*1
AnyC = playerTable.Position.str.contains('C')*1
AnyPos = [1]*500
ConVec_Position = pd.DataFrame({'AnyPG': AnyPG, 'AnySG': AnySG, 'AnySF': AnySF, 'AnyPF': AnyPF, 'AnyC': AnyC, 'AnyPos': AnyPos})

# Make salary constraint vectors
ConVec_Salary = playerTable[['Salary']]

# Get Objective (projected points)
ProjectedPoints = playerTable['ProjPts'].tolist()

# Make constraint values and directions
ConValGT = [1,1,1,1,1, 8] # Must have at least one of each position, and 8 total
ConValLT = [3,3,3,3,2, 8, 50000] # Can't have more than 3 of each position, 8 total, and 50000 salary
ConDirGT = [1] * (len(ConVec_Position.columns)) # Make >= directions for each position constraint
ConDirLT = [-1] * (len(ConVec_Position.columns) + 1) # Make <= directions for each position constraint, plus 1 for Salary

# Compile constraint vectors, values, and directions
ConVec_Final = pd.concat([ConVec_Position, ConVec_Position, ConVec_Salary], axis=1)
ConVal_Final = ConValGT + ConValLT
ConDir_Final = ConDirGT + ConDirLT

# Force all decisions to be binary
vLB = [0] * len(ProjectedPoints) # lower bound 0
vUB = [1] * len(ProjectedPoints) # upper bound 1
xint = [i+1 for i in range(len(ProjectedPoints))] # all decisions are integers (aka all are binary)

solveTimefix = time.time()
[obj, sol, duals] = lp_solve(ProjectedPoints, 
                             ConVec_Final.T.values.tolist(), 
                             ConVal_Final, 
                             ConDir_Final, 
                             vLB, 
                             vUB, 
                             xint)

solveTimeCurrRun = (time.time() - solveTimefix)
print("\n--- It took %s seconds to get all LUs ---
" % round(solveTimeCurrRun, 2))

1. 更换求解器包：使用pulp调用同源lp_solve #

R版用的lpSolve和Python的pulp可以调用同一底层lp_solve库，消除不同封装层的性能差异。优化后代码：

import pandas as pd
import numpy as np
import pulp
import time

# 生成数据（保留原逻辑）
ID = range(500)
Positions = ['PG', 'SG', 'SF', 'PF', 'C', 'PG/SG', 'SG/SF', 'PF/C']
Position = np.random.choice(Positions,500)
Salary = np.random.choice(range(3000,12000,100),500)
playerTable = pd.DataFrame({'ID': ID, 'Position': Position, 'Salary': Salary})
playerTable['ProjPts'] = playerTable['Salary']/1000*5

start_time = time.time()
# 初始化最大化模型
model = pulp.LpProblem("FantasyBasketball", pulp.LpMaximize)

# 创建二进制决策变量：每个球员选/不选
players = pulp.LpVariable.dicts("Player", playerTable.index, lowBound=0, upBound=1, cat='Binary')

# 目标函数：最大化总预测得分
model += pulp.lpSum([playerTable.loc[i, 'ProjPts'] * players[i] for i in playerTable.index])

# 位置约束：至少1人
model += pulp.lpSum([players[i] for i in playerTable.index if 'PG' in playerTable.loc[i, 'Position']]) >= 1
model += pulp.lpSum([players[i] for i in playerTable.index if 'SG' in playerTable.loc[i, 'Position']]) >= 1
model += pulp.lpSum([players[i] for i in playerTable.index if 'SF' in playerTable.loc[i, 'Position']]) >= 1
model += pulp.lpSum([players[i] for i in playerTable.index if 'PF' in playerTable.loc[i, 'Position']]) >= 1
model += pulp.lpSum([players[i] for i in playerTable.index if 'C' in playerTable.loc[i, 'Position']]) >= 1

# 位置约束：最多人数
model += pulp.lpSum([players[i] for i in playerTable.index if 'PG' in playerTable.loc[i, 'Position']]) <= 3
model += pulp.lpSum([players[i] for i in playerTable.index if 'SG' in playerTable.loc[i, 'Position']]) <= 3
model += pulp.lpSum([players[i] for i in playerTable.index if 'SF' in playerTable.loc[i, 'Position']]) <= 3
model += pulp.lpSum([players[i] for i in playerTable.index if 'PF' in playerTable.loc[i, 'Position']]) <= 3
model += pulp.lpSum([players[i] for i in playerTable.index if 'C' in playerTable.loc[i, 'Position']]) <= 2

# 总人数约束：恰好8人
model += pulp.lpSum([players[i] for i in playerTable.index]) == 8

# 薪资约束：总薪资不超50000
model += pulp.lpSum([playerTable.loc[i, 'Salary'] * players[i] for i in playerTable.index]) <= 50000

# 关闭日志输出，调用lp_solve求解
model.solve(pulp.LpSolverDefault(msg=0))

solve_time = round(time.time() - start_time, 4)
print(f"\n--- It took {solve_time} seconds to optimize lineups ---")

2. 消除冗余数据结构转换 #

原代码中频繁在DataFrame、列表、矩阵间转换带来额外开销，优化思路：

• 用numpy数组预计算位置掩码，避免重复调用str.contains
• 直接构建约束矩阵，减少不必要的拼接和转置

# 预计算位置掩码（numpy数组，一次计算复用）
pg_mask = np.char.find(playerTable['Position'].values.astype(str), 'PG') != -1
sg_mask = np.char.find(playerTable['Position'].values.astype(str), 'SG') != -1
sf_mask = np.char.find(playerTable['Position'].values.astype(str), 'SF') != -1
pf_mask = np.char.find(playerTable['Position'].values.astype(str), 'PF') != -1
c_mask = np.char.find(playerTable['Position'].values.astype(str), 'C') != -1
all_mask = np.ones(500, dtype=int)

# 直接构建约束矩阵（numpy数组，仅转置一次）
constraint_matrix = np.vstack([
    pg_mask, sg_mask, sf_mask, pf_mask, c_mask, all_mask,
    pg_mask, sg_mask, sf_mask, pf_mask, c_mask, all_mask, playerTable['Salary'].values
]).astype(int)

3. 配置求解器加速参数 #

关闭冗余输出、设置启发式优先级，进一步压缩求解时间：

• 在pulp中通过msg=0关闭日志
• 对于lp_solve，可设置timeout、presolve等参数提前终止非必要计算

4. 预缓存静态结构（针对数千次运行） #

如果要重复运行数千次，提前生成并缓存位置掩码、约束矩阵模板等静态结构，避免每次运行重复计算，可大幅降低总耗时。

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