如何加速大量Pulp优化函数调用？性能优化问询

问题背景 #

现有一个基于PuLP实现的库存买卖收益最大化函数：

• 基础场景：处理2组20维1D数组，返回20维数组，运行正常
• 进阶场景：处理20×100的2D数组，通过map循环调用100次，返回20×100数组
• 当前痛点：处理N×20×100的3D数组（对应N天、100个场景）时，嵌套循环调用性能极差，需要优化执行效率，预期返回N×20×2×100维度的结果数组。

附核心原始函数代码：

# ----------- IMPORTS --------------------------------------------------
import numpy as np
from pulp import *
# ----------- END OF IMPORTS ---------------------------------------------


def revenue_max(buy_prices, sell_prices, num_months, capacity, max_buy, max_sell, start_inventory, reqd_end_inventory):
    """
    Optimizes buy and sell quantities for next budget period
    Parameters:
    buy_prices: num_months x 1 array of floats; 
    sell_prices: num_months x 1 array of floats; 
    capacity: decimal integer; inventory capacity
    max_buy: decimal integer; max for each month
    max_sell: decimal integer; max for each month
    start_inventory: anticipated inventory at start of budget period
    reqd_end_inventory: requirement on budget period end inventory
    Constraints:
    inventory: nx1 floats; inventory >= 0; inventory <= capacity
    buy: optimal buy for each month; buy >= 0; buy <= max_buy
    sell: optimal sell  for each month; sell >= 0; sell <= max_sell
    buy and sell same month not allowed
    Returns:
    all_data: num_months x 3 array, contains buy/sell quantities, cash flows
    """

    M = 1000000     # BIG M
    all_data = np.zeros((num_months, 3))

    # ###################################
    period_start_inventory = np.zeros(num_months)
    period_end_inventory = np.zeros(num_months)
    daily_end_inv = np.zeros(shape=(num_months, 31))

    # Creating the lists to be exported to MS Excel file through pandas dataframe
    buyValue = list()
    sellValue = list()
    iValue = list()
    gValue = list()

    budget_horizon = list(range(num_months))
    # #########################Defining the variables########################
    buy = LpVariable.matrix("buy", budget_horizon, 0, None,
                            LpInteger)  # Variable-1 to 24.Injection quantity during time t
    sell = LpVariable.matrix("sell", budget_horizon, 0, None,
                             LpInteger)  # Variable-25 to 48.Withdrawal quantity during time t
    inv = LpVariable.matrix("inv", budget_horizon, 0, None,
                            LpInteger)  # Variable-49 to 72.Inventory in hand at the beginning of time t
    z = LpVariable.matrix("z", budget_horizon, 0, 1,
                          LpBinary)  # Variable-73 to 96.Inventory in hand at the beginning of time t
    # ##########################Objective function (Max. expected profit)############################
    prob = LpProblem(name="Buy_Sell_Schedule", sense=LpMaximize)
    prob += (lpSum([-buy_prices[t] * buy[t] +
                    sell_prices[t] * sell[t]
                    for t in
                    budget_horizon]))
    # ############################Defining the constraints#########################################
    for t in list(range(num_months)):
        prob += sell[t] <= M * z[t]
        prob += buy[t] <= M * (1 - z[t])
        prob += inv[t] <= capacity

        if t > 0:
            b = t - 1
            # buy
            prob += (buy[t] <= capacity - inv[b])
            # sell
            prob += (sell[t] <= inv[b])
            # inventory
            prob += inv[t] == inv[b] + buy[t] - sell[t]  # ORIGINAL

    # Period-0 buy constraint
    prob += (buy[0] <= capacity - start_inventory)

    # Period-0 sell constraint
    prob += (sell[0] <= start_inventory)

    # Period-0 inventory constraint
    prob += (inv[0] == buy[0] - sell[0] + start_inventory)

    # last budget period month constraints
    prob += inv[num_months - 1] == inv[num_months - 2] + buy[num_months - 1] - sell[
        num_months - 1]
    prob += inv[num_months - 1] >= reqd_end_inventory
    prob += inv[num_months - 1] <= reqd_end_inventory

    # ###############################Solving with default cbc solver###################################
    prob.solve(PULP_CBC_CMD(msg=False))
    # prob.solve(GLPK_CMD(msg=False))
    # ################################################t##################
    # Extracting the required variables from the model output
    for t in list(range(num_months)):  # budget_horizon:
        buyValue.append(str(buy[t].varValue))  # Injection quantities
        sellValue.append(str(sell[t].varValue))  # withdrawal quantities
        if t == 0:
            inv[t] = (buy[t].varValue - sell[t].varValue + start_inventory)
        else:
            inv[t] = inv[t - 1] + (buy[t].varValue - sell[t].varValue)
        iValue.append(inv)
        g = (-buy[t].varValue * buy_prices[t] + sell[t].varValue * sell_prices[t])  # Monthly cashflow
        gValue.append(g)
       
    buy_qty = np.array(buyValue, dtype=float)
    sell_qty = np.array(sellValue, dtype=float)
    schedule_price = np.where(buy_qty > 0, buy_prices,
                           np.where(sell_qty > 0, sell_prices, 0))

    # --- start and end inventory for each month
    period_start_inventory[0] = start_inventory
    period_end_inventory[0] = period_start_inventory[0] + buy_qty[0] - sell_qty[0]
    for j in range(1, num_months):
        period_start_inventory[j] = period_end_inventory[j - 1]
        period_end_inventory[j] = period_start_inventory[j] + buy_qty[j] - sell_qty[j]

    buy_sell_qty = np.where((buy_qty > 0), -buy_qty, np.where(sell_qty > 0, sell_qty, 0))

    # CASH_FLOWS
    # --- multiply volumes by schedule prices to get cash flows
    monthly_cash_flows_schedule = buy_sell_qty * schedule_price

    # COLLECT ALL RETURN DATA ---- dimension = num_months x 1
    all_data[:, 0] = buy_sell_qty  
    all_data[:, 1] = monthly_cash_flows_schedule

    return all_data

优化方案 #

1. 并行化处理（最快见效） #

每个场景（N×20×100中的第三个维度）的求解完全独立，适合用多进程并行替代串行循环/单线程map，充分利用CPU多核资源。

实现代码：

import numpy as np
from pulp import *
from concurrent.futures import ProcessPoolExecutor

# 先优化单个函数的冗余操作
def optimized_revenue_max(args):
    buy_prices, sell_prices, num_months, capacity, max_buy, max_sell, start_inventory, reqd_end_inventory = args
    M = 1000000
    all_data = np.zeros((num_months, 2))  # 只保留需要的两列，去掉冗余的第三列

    budget_horizon = range(num_months)
    # 定义变量
    buy = LpVariable.matrix("buy", budget_horizon, lowBound=0, cat=LpInteger)
    sell = LpVariable.matrix("sell", budget_horizon, lowBound=0, cat=LpInteger)
    inv = LpVariable.matrix("inv", budget_horizon, lowBound=0, upBound=capacity, cat=LpInteger)
    z = LpVariable.matrix("z", budget_horizon, lowBound=0, upBound=1, cat=LpBinary)

    # 目标函数
    prob = LpProblem("Buy_Sell_Schedule", LpMaximize)
    prob += lpSum(-buy_prices[t] * buy[t] + sell_prices[t] * sell[t] for t in budget_horizon)

    # 约束条件
    for t in budget_horizon:
        # 同一月不能同时买卖
        prob += sell[t] <= M * z[t]
        prob += buy[t] <= M * (1 - z[t])

        if t > 0:
            prev_inv = inv[t-1]
            prob += buy[t] <= capacity - prev_inv
            prob += sell[t] <= prev_inv
            prob += inv[t] == prev_inv + buy[t] - sell[t]

    # 初始月约束
    prob += buy[0] <= capacity - start_inventory
    prob += sell[0] <= start_inventory
    prob += inv[0] == start_inventory + buy[0] - sell[0]

    # 期末库存约束
    prob += inv[-1] == reqd_end_inventory

    # 求解
    prob.solve(PULP_CBC_CMD(msg=False))

    # 提取结果（简化操作，避免字符串转换）
    buy_qty = np.array([v.varValue for v in buy], dtype=float)
    sell_qty = np.array([v.varValue for v in sell], dtype=float)

    # 计算买卖量和现金流
    buy_sell_qty = np.where(buy_qty > 0, -buy_qty, np.where(sell_qty > 0, sell_qty, 0))
    schedule_price = np.where(buy_qty > 0, buy_prices, np.where(sell_qty > 0, sell_prices, 0))
    monthly_cash_flows = buy_sell_qty * schedule_price

    all_data[:, 0] = buy_sell_qty
    all_data[:, 1] = monthly_cash_flows
    return all_data

# 处理3D数组的函数
def process_3d_input(buy_prices_3d, sell_prices_3d, capacity, max_buy, max_sell, start_inventory, reqd_end_inventory):
    """
    buy_prices_3d: N×20×100的3D数组，对应N天、20个月、100个场景的买入价格
    sell_prices_3d: N×20×100的3D数组，对应N天、20个月、100个场景的卖出价格
    返回: N×20×2×100的结果数组
    """
    N_days, num_months, N_scenarios = buy_prices_3d.shape
    result = np.zeros((N_days, num_months, 2, N_scenarios))

    # 遍历每一天，并行处理当天的所有场景
    for day_idx in range(N_days):
        # 准备当天所有场景的参数列表
        args_list = []
        for scenario_idx in range(N_scenarios):
            buy_p = buy_prices_3d[day_idx, :, scenario_idx]
            sell_p = sell_prices_3d[day_idx, :, scenario_idx]
            args = (buy_p, sell_p, num_months, capacity, max_buy, max_sell, start_inventory, reqd_end_inventory)
            args_list.append(args)
        
        # 并行求解当天的所有场景
        with ProcessPoolExecutor() as executor:
            day_results = list(executor.map(optimized_revenue_max, args_list))
        
        # 将结果填入3D数组
        for scenario_idx, res in enumerate(day_results):
            result[day_idx, :, :, scenario_idx] = res
    
    return result

2. 单个模型的细节优化 #

原函数存在不少冗余操作，优化后可减少单个模型的构建和求解时间：

• 变量定义时直接设置lowBound/upBound，减少约束数量（比如库存上限可直接在inv变量定义中设置）
• 去掉无用的变量（如daily_end_inv、iValue、gValue等未用到的列表）
• 结果提取时直接获取变量值，避免字符串转换再转float
• 简化期末库存约束，合并>=和<=为==

3. 批量建模（可选，复杂度较高） #

如果并行化仍不够快，可以尝试将多个场景的模型合并成一个大模型求解：

• 为每个场景的变量添加前缀（如buy_s0_t0）
• 复制约束到每个场景
• 目标函数为所有场景的收益之和
• 求解一次大模型后拆分结果

这种方法适合场景数量不多的情况，变量规模过大时求解器效率会下降。

注意事项 #

• 多进程并行时，Windows系统下需将主逻辑放在if __name__ == '__main__':块中，避免子进程重复初始化
• 可尝试更换商业求解器（如Gurobi、CPLEX），其在大规模问题上性能远优于默认的CBC求解器
• 若场景参数存在重复，可缓存已求解的结果，避免重复计算

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