多彩编程 多彩编程MZPH · CODE BLOG
ARTICLE DETAIL

文章详情

深耕前端与后端开发技术的一线实战笔记与踩坑复盘。

国赛C题供应链优化:基于MILP的随机需求订购与运输联合决策模型

国赛C题供应链优化:基于MILP的随机需求订购与运输联合决策模型 1. 项目概述从“问题二”看国赛C题的实战精髓每年九月的全国大学生数学建模竞赛国赛对很多理工科学生来说都是一场脑力与耐力的双重考验。2021年的C题以其贴近实际、数据复杂、模型综合的特点给参赛队伍留下了深刻印象。题目聚焦于“生产企业原材料的订购与运输”这是一个典型的供应链优化问题而其中的“问题二”更是整个赛题承上启下的关键枢纽。它不再仅仅是简单的计算而是要求我们基于第一问建立的订购方案进一步考虑运输成本、库存策略与随机需求的动态博弈构建一个完整的、成本最优的订购与运输决策模型。简单来说如果第一问是让你学会“买多少”那么问题二就是考验你如何“聪明地买”和“高效地运”。它模拟了一个现实生产经理面临的真实困境供应商的批发价有折扣门槛运输费用按车次计费且存在“拼车”优惠未来需求存在不确定性仓库容量还有限制。你需要在这些错综复杂的约束条件下找到一个成本包含订购费、运输费、库存费和可能的缺货损失最低的运营策略。这不仅仅是一个数学问题更是一个需要融合运筹学、统计学和计算机仿真的系统工程问题。对于参赛者而言攻克问题二意味着掌握了解决一类经典库存与运输联合优化问题的核心方法论。无论是使用MATLAB进行矩阵运算和规划求解还是利用Python的Pandas、NumPy进行数据处理配合Scikit-learn或Statsmodels进行预测亦或是深入使用Gurobi、PuLP等优化求解器这道题都提供了一个绝佳的练兵场。接下来我将以一名多次参与竞赛指导的“老队员”视角拆解问题二的完整解决思路、核心模型构建、算法实现细节以及那些在论文里不会写的“踩坑”实录。2. 核心思路拆解如何构建订购与运输的联合优化模型面对问题二最忌讳的就是一头扎进代码里。清晰的思路分层是成功的一半。我的策略通常分为四个层次理解约束、定义决策、构建目标、选择算法。2.1 问题约束的深度解读题目给出的约束就是游戏的规则必须逐字逐句吃透原材料订购约束供应商的销售策略是阶梯价格。例如当周购买量低于某个值A时单价为P1超过A但低于B时超过A的部分单价为P2P2 P1超过B时超过B的部分单价为P3P3 P2。这直接决定了我们的订购成本函数是一个分段线性函数。运输约束这是本题最大的难点和特色之一。运输公司按“车”计费每车有固定容量如6立方米。但运费不是简单的“车数×单价”而是有“拼车”优惠当本周运输量大于等于3车时运费打9折。这意味着运输成本函数也是一个分段函数且决策变量运输量与成本之间的关系是非线性的、离散的。库存动态与容量约束仓库有最大容量限制。每周的库存状态是动态变化的本周库存 上周库存 本周到货量 - 本周消耗量。消耗量由生产需求决定而需求是随机的题目附件会给出历史数据或分布。库存不能为负但允许缺货缺货可能产生惩罚成本需根据题目具体说明确定也不能超过仓库上限。需求不确定性未来每周的需求是随机变量。我们需要处理这种不确定性常见的思路有两种一是将随机需求用其期望值代替转化为确定性模型二是采用随机规划或模拟仿真的方法如使用情景分析法Scenario Analysis或蒙特卡洛模拟。2.2 决策变量与目标函数的形式化基于约束我们需要定义清晰的决策变量。对于未来N周例如24周的计划期我们需要决定x_t: 第t周的原材料订购量立方米。实际上由于运输按车计更精细的变量可以是y_t: 第t周租用的运输车辆数量整数。那么运输量就是y_t * 单车容量。 但订购量x_t不一定等于当周到货量这里可能涉及运输决策。为了简化许多优秀论文假设订购后立即运输即当周到货那么x_t同时也是运输量。我们采用这个简化假设进行阐述。因此我们的核心决策变量是一个序列{x_1, x_2, ..., x_N}其中每个x_t是连续变量订购量但与之相关的运输成本又依赖于离散的车辆数。目标函数是最小化总期望成本它由四部分组成订购成本 C_order(x_t)由阶梯价格规则决定的分段线性函数。运输成本 C_transport(x_t)由车辆数和折扣规则决定的分段函数。车数 ceil(x_t / 单车容量)其中ceil是向上取整函数。库存持有成本 C_holding(I_t)I_t * 单位库存持有费其中I_t是第t周末的库存水平。缺货惩罚成本 C_shortage(S_t)max(0, -I_t) * 单位缺货惩罚费如果允许缺货且有惩罚。总目标Minimize E[ Σ_{t1}^{N} (C_order(x_t) C_transport(x_t) C_holding(I_t) C_shortage(S_t)) ]其中期望E是针对未来随机需求序列的。I_t由库存动态方程I_t I_{t-1} x_t - d_t决定d_t是随机需求。2.3 模型类型的判断与算法选型这是一个典型的动态的、随机的、带有非线性分段成本函数和整数约束车辆数的库存控制问题。直接求解析解几乎不可能必须借助数值优化方法。主流思路有两种思路一随机动态规划SDP或模型预测控制MPC。将多期问题分解为一系列单期问题。在每一期初根据当前库存I_{t-1}和未来需求的预测分布求解一个单期或有限期如滚动3期的优化问题得到当前的最优订购量x_t。执行x_t后需求实现库存更新进入下一期。这种方法能较好地处理不确定性但计算量较大需要对需求分布进行估计或采用情景树。思路二转化为确定型混合整数线性规划MILP进行近似求解。这是比赛中更实用、更易实现的方法。其核心技巧在于线性化线性化阶梯价格引入0-1辅助变量z_{t1}, z_{t2}, z_{t3}表示落在哪个价格区间并引入连续变量表示各区间的采购量通过大M法构造线性约束。线性化运输折扣引入0-1辅助变量w_t当运输车数≥3时w_t1否则为0。运输成本表示为[车数 * 单价 * (1 - 折扣率 * w_t)]同样用大M法将车数与w_t的关系转化为线性约束。处理随机需求采用情景分析法。根据历史数据用时间序列模型如ARIMA或分布拟合生成M个可能的需求情景每个情景是一个长度为N的需求序列。目标函数变为最小化所有情景下的平均总成本。这样随机规划就转化为了一个大规模但确定性的MILP问题。对于国赛有限的时间思路二情景分析MILP的可行性最高。它思路直观有成熟的求解器如Gurobi, CPLEX支持论文中也易于阐述。我们后续的代码实现也将围绕此思路展开。3. 关键步骤实现与代码解析Python为例这里我以Python生态为例展示如何将上述思路落地。我们会用到pandas处理数据numpy进行计算scikit-learn或statsmodels进行需求预测以生成情景最后用pulp或gurobipy来构建和求解MILP模型。3.1 数据预处理与需求情景生成假设我们已从附件中读取了历史需求数据demand_history。import pandas as pd import numpy as np from statsmodels.tsa.arima.model import ARIMA import itertools # 1. 读取与预处理数据 # df pd.read_excel(附件.xlsx) # demand_history df[历史需求].values # 这里我们用模拟数据 np.random.seed(2021) demand_history np.random.normal(loc100, scale20, size100) # 100周历史数据 # 2. 利用ARIMA模型拟合并生成未来N周的多组情景 def generate_demand_scenarios(history, n_periods24, n_scenarios50, order(1,1,1)): 生成未来需求的情景。 history: 历史需求序列 n_periods: 需要预测的未来期数 n_scenarios: 生成的情景数量 order: ARIMA模型阶数可通过AIC/BIC准则选择 model ARIMA(history, orderorder) model_fit model.fit() scenarios [] for _ in range(n_scenarios): # 使用拟合的模型进行模拟生成一条未来路径 forecast model_fit.simulate(n_periods, anchorend) # 确保需求非负 forecast np.maximum(forecast, 0) scenarios.append(forecast) return np.array(scenarios) # 形状: (n_scenarios, n_periods) N 24 # 规划期24周 S 50 # 生成50个需求情景 demand_scenarios generate_demand_scenarios(demand_history, N, S, order(1,1,0)) print(f生成了{S}个未来{N}周的需求情景形状{demand_scenarios.shape})注意ARIMA模型的选择和参数调优需要时间。比赛中如果数据特征明显如季节性可尝试SARIMA或更简单的移动平均、指数平滑。如果时间紧迫一个实用的“土办法”是对历史数据做Bootstrap重采样随机组合出未来情景这也能捕捉历史波动性。3.2 构建混合整数线性规划MILP模型我们使用pulp库它是一个免费的线性规划建模接口。import pulp as pl # 定义问题参数示例值需根据题目更改 C 6.0 # 单车容量立方米 unit_trans_cost 1000.0 # 单车运输费用元/车 discount_rate 0.1 # 运输量3车时的折扣率 (9折) discount_threshold 3 # 享受折扣的最低车数 # 阶梯价格参数 price_breaks [0, 200, 500, float(inf)] # 数量分段点 unit_prices [0, 1000, 950, 900] # 各分段单价第一个为0占位 # 库存参数 init_inv 50 # 期初库存立方米 max_inv 1000 # 最大库存容量 holding_cost 50 # 单位库存持有成本元/立方米/周 shortage_cost 200 # 单位缺货惩罚成本元/立方米/周若不允许缺货则设为极大值 # 创建问题 prob pl.LpProblem(Material_Ordering_Transport_Problem, pl.LpMinimize) # 定义决策变量 # 连续变量每周的订购量 x {t: pl.LpVariable(fx_{t}, lowBound0, catContinuous) for t in range(1, N1)} # 整数变量每周的运输车数 y {t: pl.LpVariable(fy_{t}, lowBound0, catInteger) for t in range(1, N1)} # 0-1变量是否享受运输折扣 w {t: pl.LpVariable(fw_{t}, catBinary) for t in range(1, N1)} # 阶梯价格相关的辅助变量以第一周t1为例实际需要循环所有t和所有分段k # 这里展示处理一个分段区间例如第2段200-500的通用方法实际需要为每个t定义多组变量 # 假设我们定义三个分段区间 (0-200], (200-500], (500-∞) # 引入0-1变量z_{t,k}表示第t周的订购量是否落在第k个分段 z {} for t in range(1, N1): for k in range(1, len(price_breaks)): # k1,2,3 代表三个分段 z[(t, k)] pl.LpVariable(fz_{t}_{k}, catBinary) # 引入连续变量a_{t,k}表示第t周在第k个分段的订购量 a {} for t in range(1, N1): for k in range(1, len(price_breaks)): a[(t, k)] pl.LpVariable(fa_{t}_{k}, lowBound0, catContinuous) # 库存变量每个情景s下每周的库存水平 I {(s, t): pl.LpVariable(fI_{s}_{t}, lowBound-float(inf), catContinuous) for s in range(S) for t in range(0, N1)} # 缺货变量每个情景s下每周的缺货量 short {(s, t): pl.LpVariable(fshort_{s}_{t}, lowBound0, catContinuous) for s in range(S) for t in range(1, N1)} # 持货变量每个情景s下每周的库存持有量 hold {(s, t): pl.LpVariable(fhold_{s}_{t}, lowBound0, catContinuous) for s in range(S) for t in range(1, N1)} # 添加约束 M 10000 # 一个大数用于大M法线性化 # 约束1: 订购量与分段变量的关系 for t in range(1, N1): # 订购总量等于各分段量之和 prob x[t] pl.lpSum(a[(t, k)] for k in range(1, len(price_breaks))) # 每个分段量a_{t,k}的范围约束大M法 for k in range(1, len(price_breaks)): lower price_breaks[k-1] upper price_breaks[k] prob a[(t, k)] lower * z[(t, k)] prob a[(t, k)] upper * z[(t, k)] # 确保只有一个分段被激活 prob pl.lpSum(z[(t, k)] for k in range(1, len(price_breaks))) 1 # 约束2: 运输车数与订购量的关系 (y_t ceil(x_t / C)) 的线性近似 # 我们将其松弛为 x_t y_t * C 并且 x_t (y_t - 1) * C # 由于y是整数这等价于 y_t ceil(x_t / C)。在最小化成本运输成本随y_t增加的目标下求解器会自动使y_t取满足条件的最小整数。 for t in range(1, N1): prob x[t] y[t] * C prob x[t] (y[t] - 1) * C - M * (1 - w[t]) # 这是一个近似严格线性化需要引入额外变量此处简化处理。更严谨的做法是引入另一个0-1变量处理“大于”关系。 # 约束3: 运输折扣逻辑线性化 for t in range(1, N1): # 如果 y_t discount_threshold, 则 w_t 可以被强制为1享受折扣 prob y[t] discount_threshold * w[t] # 如果 y_t discount_threshold, 则 w_t 可以被强制为0 prob y[t] discount_threshold - 1 M * w[t] # 注意上述约束是“或”的关系标准大M法需要更复杂的构造。一个更简洁但近似的做法是 # prob y[t] discount_threshold * w[t] # prob y[t] M * w[t] (discount_threshold - 1) * (1 - w[t]) # 这个约束更准确 # 约束4: 库存动态平衡对每个情景s for s in range(S): # 初始库存 prob I[(s, 0)] init_inv for t in range(1, N1): # 库存平衡方程: I_{s,t} I_{s,t-1} x_t - demand_scenarios[s, t-1] prob I[(s, t)] I[(s, t-1)] x[t] - demand_scenarios[s, t-1] # 库存容量约束 prob I[(s, t)] max_inv # 将库存水平分解为持货和缺货部分: I_{s,t} hold_{s,t} - short_{s,t} prob I[(s, t)] hold[(s, t)] - short[(s, t)] # 持货和缺货非负 prob hold[(s, t)] 0 prob short[(s, t)] 0 # 定义目标函数最小化所有情景下的平均周成本 order_trans_cost pl.lpSum([ pl.lpSum(unit_prices[k] * a[(t, k)] for k in range(1, len(price_breaks))) # 订购成本 unit_trans_cost * y[t] * (1 - discount_rate * w[t]) # 运输成本 for t in range(1, N1) ]) inventory_cost (1/S) * pl.lpSum([ # 求所有情景的平均 holding_cost * hold[(s, t)] shortage_cost * short[(s, t)] for s in range(S) for t in range(1, N1) ]) prob order_trans_cost inventory_cost # 求解问题 solver pl.PULP_CBC_CMD(msgFalse, timeLimit300) # 使用CBC求解器静默模式限制5分钟 # 如果有Gurobi可以替换为 solver pl.GUROBI_CMD(timeLimit300) prob.solve(solver) # 输出求解状态和结果 print(f求解状态: {pl.LpStatus[prob.status]}) print(f最优总期望成本: {pl.value(prob.objective):.2f} 元) if prob.status pl.LpOptimal: print(\n前4周的最优订购方案示例:) for t in range(1, 5): print(f 第{t}周: 订购量{pl.value(x[t]):.2f}, 车数{pl.value(y[t])}, 是否折扣{是 if pl.value(w[t])0.5 else 否})这段代码构建了一个完整的、考虑随机需求的订购与运输联合优化模型。它包含了阶梯价格、运输折扣、库存动态和容量约束的线性化处理。4. 模型求解的优化技巧与注意事项直接运行上述模型在情景数S50、周期N24时变量和约束数量已经非常庞大数千个变量和约束对免费求解器CBC是个挑战。在实际比赛中需要一些技巧来提升求解效率和方案质量。4.1 求解效率提升策略减少情景数量50个情景可能过多。可以通过聚类算法如K-Means对生成的大量原始情景进行聚类选取5-10个具有代表性的“典型情景”及其概率权重进行计算。这能大幅降低模型规模。from sklearn.cluster import KMeans # 假设demand_scenarios形状为 (500, 24)我们先生成500条路径 kmeans KMeans(n_clusters10, random_state0).fit(demand_scenarios) cluster_centers kmeans.cluster_centers_ # 10个典型情景 cluster_weights np.bincount(kmeans.labels_) / len(kmeans.labels_) # 对应概率 # 然后用这10个加权情景代替原来的50个简化模型运输折扣的精确线性化前面代码中的折扣约束是近似。更精确的线性化需要引入额外的连续变量和约束但这会增加复杂度。如果时间紧可以放弃精确线性化采用两阶段法先求解不考虑折扣的模型得到订购量序列x_t再根据x_t计算最优的运输车辆安排y_t和折扣决策w_t。虽然可能不是全局最优但结果合理且计算快。松弛整数变量可以先求解y_t为连续变量的松弛问题得到解后对y_t向上取整作为初始解提供给MILP求解器能加速求解。利用求解器特性如果使用Gurobi可以设置MIPGap允许的间隙为一个较小的值如0.01这样求解器在找到可行解与最优解差距在1%以内时就会停止节省时间。在Pulp中可以通过参数设置。4.2 方案分析与灵敏度测试得到最优解后不能只报一个数字了事。你需要分析方案的特征这往往是论文的加分项。策略分析观察最优的x_t序列。它是否呈现某种规律例如是否在需求预测高的前期提前备货追逐策略还是保持相对平稳的订购量均衡策略结合库存水平I_t的变化图进行分析。成本构成分析计算总成本中订购费、运输费、库存持有费和缺货费的占比。这能揭示成本控制的重点。例如如果运输费占比极高说明“拼车”策略运用得不够可能需要调整模型或参数。灵敏度分析改变关键参数观察方案和总成本的变化。这是体现模型稳健性的关键。需求波动性增大生成情景时的方差看最优成本增加多少方案是否变化剧烈。价格参数如果供应商提高折扣门槛如从200提高到250总成本会增加多少库存成本如果仓库租金上涨holding_cost增加最优策略是否会倾向于更频繁的小批量订购Just-in-Time缺货惩罚如果缺货导致停产损失巨大shortage_cost极高策略是否会变得非常保守维持很高的安全库存进行灵敏度分析后你可以给出管理启示例如“建议企业与运输公司谈判争取将享受折扣的车数门槛从3车降低到2车预计可降低总成本约X%”。5. 常见“踩坑点”与实战心得回顾多次竞赛和指导经历队伍在解决这类问题时最容易在以下几个地方“翻车”。对“运输折扣”建模错误这是最高频的错误。很多人直接用if y_t 3: cost ... else: cost ...的逻辑写进目标函数但这在MILP中是非线性的。必须使用0-1变量进行线性化。另一个常见错误是忽略了“车数”是整数直接用x_t / C计算成本。心得遇到“如果...则...”的成本或约束第一时间想到引入0-1辅助变量和大M法线性化。大M的取值要足够大以起到约束作用但又不能太大否则会影响求解精度一般取一个比该变量可能最大值稍大的数即可。库存动态方程符号混乱I_t I_{t-1} x_t - d_t中的d_t是“消耗”还是“需求”如果d_t是生产需求那么它就是出库量。务必明确每个变量的物理意义并在论文中给出清晰的定义。符号混乱会导致整个模型全错。忽略初始库存和期末库存处理期初库存I_0是已知参数。对于期末库存I_N题目有时会要求不低于某个安全库存或者对期末库存本身有成本计算如期末剩余原材料按价值折旧处理。务必仔细阅读题目对首尾周的特殊说明。随机处理过于简单或过于复杂有的队伍直接用历史平均值作为未来各周需求这完全忽略了不确定性模型结果过于乐观。有的队伍则执着于构建极其复杂的随机过程模型消耗大量时间却收效甚微。心得“中庸之道”最实用。情景分析法是平衡准确性与复杂度的最佳选择。生成情景时如果历史数据少可以假设需求服从正态分布需检验用均值和方差生成随机数如果数据有一定趋势或季节性用简单的指数平滑或ARIMA(1,1,0)足矣。重点是把更多时间花在优化模型的构建和求解上。代码调试耗时过长在编写复杂的MILP模型时很容易出现约束逻辑错误导致模型不可行或无界。调试技巧从简到繁逐步验证。先构建一个只有1个情景、2个周期的最简化模型手动计算一个预期结果看求解器输出是否匹配。然后逐步增加周期数、情景数、引入分段价格、引入运输折扣等复杂要素。每增加一个功能都检查模型状态和结果是否合理。另外将模型输出prob.writeLP(model.lp)写入文件可以直观检查所有变量和约束。论文表述与模型脱节论文中描述的模型和实际代码实现的模型不一致。特别是线性化部分论文里要用数学公式清晰地写出如何引入0-1变量和大M而不仅仅是文字描述。写作建议在论文的“模型建立”部分先用文字和公式定义核心决策变量、目标函数和约束。然后专门用一个小节如“5.3 非线性约束的线性化处理”来详细阐述如何将阶梯价格和运输折扣转化为线性形式。附上转化后的线性约束公式这能极大提升论文的专业性和可信度。最后记住数学建模竞赛评价的是“模型、结果、论文”三位一体。一个求解速度很快、结果看似不错的模型如果论文没有清晰表达出建模思想、求解过程和结果分析也难获高分。务必留出足够的时间将你的思路、你的代码所实现的模型、以及你从结果中挖掘的洞察清晰、准确、有条理地呈现出来。这道“订购与运输”问题本质上是在考察你如何用数学工具刻画现实世界的复杂规则并寻找最优解的能力——这种能力远比比赛本身更有价值。
返回列表