简介2022年五一赛A题《血管机器人的订购与学习优化》的完整论文与代码资料包面向数学建模竞赛参赛者、运筹优化学习者以及需要完成类似作业的学生。资源以单个PDF文件呈现文件总数仅1个、大小约1.04MB论文正文与附录代码、数据一同收录便于直接对照复现。论文围绕医院在满足每周治疗需求的同时降低运营成本建立动态规划模型用集合划分表示容器艇和操作手的数量走向从训练费、保养费、购买费等角度构造目标函数并通过Lingo求解前8周、前104周等不同阶段的最优采购方案。后续还依次引入设备损耗率、购买优惠政策以及ARIMA(3,1,4)对第105至112周需求进行预测完整覆盖建模、约束设计、求解与预测流程。附录中的代码和数据可直接用于参考仿写或课程作业该资源已有5325人学习下载能为赛题复盘和算法实践提供有力支撑。1. 血管机器人订购与学习优化先看清题目在问什么2022年五一赛A题“血管机器人的订购与学习优化”看起来是一个供应链采购问题实际上是把“买多少”和“用多久能变聪明”两件事绑在一起做联合决策。医院要采购血管介入手术机器人供应商报价不同、初始成功率不同而且机器人随着手术例数累积成功率会上升——题目里说的“学习优化”就是在采购阶段就把后续的学习曲线收益算进去让订购方案不只是便宜而是在总成本、手术需求、成功率约束之间取一个真正可落地的均衡。这篇笔记适合两类人一是数学建模竞赛参赛者想把这道题从零搭到能跑出方案的完整代码二是做医疗设备采购或供应链优化的一线工程师这套建模思路换成自己的数据一样能用。下面我按“问题拆解 → 订购建模 → 学习建模 → 联合求解 → 避坑 → 回测验证”的顺序把整条链路讲透代码全部用Python实现。2. 订购问题建模供应商、需求与成功率约束的三层结构2.1 先定变量把“订购”翻译成数学语言拿到题目后第一件事不是急着写代码而是把决策变量写清楚。我一般习惯把变量分成三层供应商层、手术类型层、时间层。常见做法是定义一套下标供应商集合S每个供应商提供不同型号的血管机器人型号之间价格、初始成功率、单台可服务的手术类型可能不同手术类型集合T例如冠状动脉介入、脑血管介入、外周血管介入等决策变量x[i][t]表示从供应商i采购、用于手术类型t的机器人台数。这里有个容易忽略的点一台机器人可能不是只能做一种手术所以x[i][t]的粒度要由题目给的数据决定。如果题目说“某型号机器人能完成多种手术”那就要拆成x[i][t]否则直接x[i]就行。实际竞赛里题目多半给出的是“供应商—型号—价格—初始成功率—是否适配某类手术”的二维表所以按x[i][t]建模更通用。目标函数自然是总采购成本最小化from pulp import LpProblem, LpMinimize, LpVariable, LpInteger, lpSum, value import pandas as pd # 示意数据供应商、手术类型、单价、适配性 suppliers [S1, S2, S3] surgery_types [T1, T2, T3] # 单价supplier x surgery_type若不适配则价格极大用于排除 price { (S1, T1): 100, (S1, T2): 110, (S1, T3): 999, (S2, T1): 95, (S2, T2): 120, (S2, T3): 105, (S3, T1): 130, (S3, T2): 999, (S3, T3): 115, } # 需求每种手术需要机器人台数 demand {T1: 12, T2: 8, T3: 10} model LpProblem(VesselRobot_Ordering, LpMinimize) x { (i, t): LpVariable(fx_{i}_{t}, 0, None, LpInteger) for i in suppliers for t in surgery_types if price[(i, t)] 900 } model lpSum(price[(i, t)] * x[(i, t)] for (i, t) in x) for t in surgery_types: model lpSum(x[(i, t)] for (i, t) in x if t t) demand[t] model.solve() print(f最低采购成本: {value(model.objective)})这段代码的逻辑是先用字典price把不适配的组合设成高价999然后在建变量时通过price 900过滤掉直接避免采购不适配的机器人。需求约束写成而不是因为多买一台虽然增加成本但可能给学习优化留出冗余后面章节会用到。关键参数是demand和x 的上界。上面的代码里没设上界求解器可能会给出非常大的整数但总成本约束会自然压住它。如果实际运行出现内存爆炸或求解过慢一定要把上界设为demand[t] * 2之类避免变量取值范围过大——这是后文避坑章的常见问题之一。2.2 成功率约束不能拍脑袋先算整体成功率光有成本还不够题目一定会给成功率要求比如“每类手术的综合成功率不得低于95%”。这里建模时最容易翻车把供应商给的初始成功率直接当最终成功率用。实际上一台机器人的成功率随手术例数增加而上升所以约束里应该用“学习后的稳态成功率”或者“预期累计成功率”。假设我们从题目数据里能拿到每个供应商的初始成功率p0[i][t]和单台最大可服务例数capacity[i][t]那么一个可用的约束是# 稳定成功率简单处理为初始成功率 学习增量 steady_rate { (S1, T1): 0.96, (S1, T2): 0.97, (S2, T1): 0.94, (S2, T2): 0.98, (S2, T3): 0.93, (S3, T1): 0.95, (S3, T3): 0.96, } min_rate 0.95 # 最低综合成功率 for t in surgery_types: total_order lpSum(x[(i, t)] for (i, t) in x if t t) if total_order 0: # 用加权平均成功率近似 model lpSum(steady_rate[(i, t)] * x[(i, t)] for (i, t) in x if t t) min_rate * total_order这里用了“加权平均成功率”近似即总成功率 每台机器人稳态成功率 × 台数之和 / 总台数。这个约束是线性的符合竞赛常用做法。更严格的做法是把学习曲线写进去得到每台机器人的最终成功率再求和但那个往往是非线性的放在下一章处理。2.3 为什么用混合整数规划而不是贪心很多参赛队上来就用贪心按单价从低到高排序先买最便宜的直到满足需求。这在小规模数据下偶尔能蒙对但一旦加入成功率约束贪心会失效——最便宜的供应商可能成功率不够导致必须额外买高价机器人来拉高加权平均这时候贪心选择的“最便宜”组合反而比多花一点钱买中价位机器人更贵。混合整数规划MIP的价值在于它能在同一个框架里同时处理“整台采购”“加权成功率”“需求满足”三种约束求解器比人脑更擅长搜索这类组合空间。我们用pulp调CBC求解器对竞赛规模的数据供应商 ≤ 20手术类型 ≤ 10通常几秒就能给出最优解。现在我们把订购问题单独跑通这只是半个题目。另一半是“学习优化”它决定steady_rate到底怎么算以及真正约束里的成功率是不是一个常数。接下来进入核心。3. 学习优化建模用贝叶斯更新修正订购决策3.1 学习曲线不是玄学是Beta分布题目里的“学习优化”常见描述是机器人每次执行手术如果成功会积累经验后续手术成功率提高如果失败也会积累教训。在工程上我们用一个Beta分布来描述一台机器人的真实成功率。Beta分布有两个参数alpha和beta分别代表“成功例数”和“失败例数”。先验可以是供应商出厂时给的数据例如某台机器人初试成功率是0.9等价于alpha9, beta1成功9次失败1次。每次手术后若成功则alpha 1失败则beta 1。这就是贝叶斯更新。更新后的期望成功率是E[p] alpha / (alpha beta)用这个期望替代上一章的steady_rate就能把学习过程融入订购模型。学习优化要回答的问题是如果新订购的机器人刚开始成功率低但随着手术例数增加会上升那么在满足总成功率要求的前提下是否值得为了长期成本而多买几台“潜力股”这直接改变了目标函数里成本与成功率的交换关系。3.2 构造学习数据并模拟手术过程竞赛题目通常会给出一个“学习效率”参数可能叫“学习因子”“经验增值”用来描述手术例数每增加一次成功率提升多少。一个常见的参数化学习曲线是p(n) p_max - (p_max - p0) * exp(-r * n)其中p0是初始成功率p_max是理论上限r是学习率n是累计手术例数。这个公式逻辑自然刚开始提升快后面逐渐平缓逼近上限。写一个Python类来模拟单个机器人的学习过程import numpy as np class LearningRobot: def __init__(self, p0, p_max, r, alpha01.0, beta01.0): self.p0 p0 self.p_max p_max self.r r # 初始Beta分布参数可由p0换算alpha0 p0*k, beta0 (1-p0)*k self.alpha alpha0 self.beta beta0 self.n 0 def current_rate(self): 用贝叶斯期望作为当前成功率估计 return self.alpha / (self.alpha self.beta) def after_n_surgeries(self, n): 返回学习曲线预测的n次手术后的期望成功率 return self.p_max - (self.p_max - self.p0) * np.exp(-self.r * n) # 示例初始成功率0.8上限0.98学习率0.05 robot LearningRobot(p00.8, p_max0.98, r0.05) for n in [0, 10, 20, 50, 100]: print(f第{n}例后理论成功率: {robot.after_n_surgeries(n):.4f})这段代码把学习曲线和贝叶斯更新分开了after_n_surgeries用于建模前的预判current_rate用于模拟过程中根据实际手术结果做后验更新。参数p0、p_max、r从哪里来如果题目给了表格直接读入如果只给了文字描述就自己设几组合理取值做灵敏度分析。注意r越大学习越快稳态来得越早。一般医学设备学习率在0.02~0.1之间属于正常范围。3.3 把学习曲线嵌入订购约束二阶段迭代法上一章的线性约束用的是固定steady_rate现在要换成学习曲线预测值。常见做法是“二阶段迭代”第一阶段用当前估计的成功率建立MIP得到一个初始订购方案。第二阶段用LearningRobot模拟该方案下每台机器人在整个任务周期内累计手术例数重新计算每台机器人结束时的期望成功率。把新的成功率带回第一阶段的约束迭代若干次直到订购方案不再变化。这种迭代和直接用非线性约束求解相比代码简单竞赛里更稳妥。尤其当题目要求“给出订购量与学习优化后的方案”迭代法能直观说明“学习优化”体现在哪里——每次迭代成功率约束的右侧会变得更宽松或更严格因为某些机器人学习后成功率提高了可能让原本无法满足95%约束的订购组合变得可行。下面给出迭代的核心循环def solve_with_learning(x, price, demand, initial_rate, learn_params): # learn_params: dict {(i,t): {p0:, p_max:, r:}} rate initial_rate.copy() for iteration in range(10): # 第一阶段按当前rate跑MIP model build_mip(x, price, demand, rate) model.solve() if model.status ! 1: # 无解或求解失败 break # 第二阶段提取订购量模拟学习 new_rate {} for (i, t), var in x.items(): qty var.varValue if var.varValue else 0 if qty 0: continue lp learn_params[(i, t)] # 假设每台机器人分到 total_demand[t] 例手术 n_per_robot demand[t] / qty new_rate[(i, t)] LearningRobot( lp[p0], lp[p_max], lp[r] ).after_n_surgeries(n_per_robot) # 更新rate只更新被订购的项未订购的保留原值 rate.update(new_rate) print(f迭代 {iteration1}: 成本{value(model.objective):.2f}, 成功率{rate}) return rate这里的build_mip是我们上一章写的模型函数唯一的改动是把steady_rate换成传入的rate。n_per_robot的估算方法很关键假设手术需求总量固定平均分摊到每台机器人。实际上不同供应商的机器人可能被分配到的手术例数不同更精细的做法是用手术类型t的总需求除以订购台数再乘以一个“利用率”系数比如0.8。你可以把利用率系数也做成参数用来调节“学习”对成功率的贡献强度。这个迭代框架的好处是你能清楚看到成功率和成本在两次迭代间的变化方便写进论文。但要注意迭代可能不收敛。下一章我们直接用完整代码把订购、学习、迭代串起来并给出一个能跑通的示例。4. 联合求解与完整代码从数据表到订购方案的最小可复现步骤4.1 准备数据自己造一张不影响逻辑的样例表竞赛题目的数据格式各不相同这里我用一张统一的表来演示你在套自己数据时只需要改csv的列名即可。表结构包含供应商ID、手术类型ID、单价、初始成功率、成功率上限、学习率、是否适配。我把不适配的单价设成极大值避免模型选中。import pandas as pd data pd.DataFrame([ [S1, T1, 100, 0.80, 0.98, 0.05, 1], [S1, T2, 110, 0.82, 0.97, 0.04, 1], [S2, T1, 95, 0.85, 0.99, 0.03, 1], [S2, T2, 120, 0.84, 0.98, 0.06, 0], # 不适配 [S2, T3, 105, 0.78, 0.96, 0.07, 1], [S3, T1, 130, 0.90, 0.99, 0.02, 1], [S3, T3, 115, 0.88, 0.98, 0.03, 1], ], columns[supplier, surgery, price, p0, pmax, r, adapter]) demand {T1: 12, T2: 8, T3: 10} min_rate 0.95把适配性做成0/1在读入时直接过滤adapter 0。单价设置为9999也可以但过滤更干净。注意我用pulp建变量时只需要进货价、初始成功率、上限、学习率这四个字段其余靠字典索引。4.2 完整代码MIP 学习迭代我把前面散落的函数合并成一个可复制的脚本。这个脚本是“示例代码讲解”级别的跑一遍就能看到订购量变化import pandas as pd import numpy as np from pulp import LpProblem, LpMinimize, LpVariable, LpInteger, lpSum, value data pd.DataFrame([ [S1, T1, 100, 0.80, 0.98, 0.05, 1], [S1, T2, 110, 0.82, 0.97, 0.04, 1], [S2, T1, 95, 0.85, 0.99, 0.03, 1], [S2, T2, 120, 0.84, 0.98, 0.06, 0], [S2, T3, 105, 0.78, 0.96, 0.07, 1], [S3, T1, 130, 0.90, 0.99, 0.02, 1], [S3, T3, 115, 0.88, 0.98, 0.03, 1], ], columns[supplier, surgery, price, p0, pmax, r, adapter]) demand {T1: 12, T2: 8, T3: 10} min_rate 0.95 # 只保留适配记录 data data[data[adapter] 1] price {(row.supplier, row.surgery): row.price for row in data.itertuples()} p0 {(row.supplier, row.surgery): row.p0 for row in data.itertuples()} pmax {(row.supplier, row.surgery): row.pmax for row in data.itertuples()} r {(row.supplier, row.surgery): row.r for row in data.itertuples()} def build_mip(rate): model LpProblem(Ordering, LpMinimize) x { (i, t): LpVariable(fx_{i}_{t}, 0, None, LpInteger) for (i, t) in price } model lpSum(price[(i, t)] * x[(i, t)] for (i, t) in x) for t in demand: model lpSum(x[(i, t)] for (i, t) in x if t t) demand[t] for t in demand: total lpSum(x[(i, t)] for (i, t) in x if t t) model lpSum(rate[(i, t)] * x[(i, t)] for (i, t) in x if t t) min_rate * total return model, x def simulate_learning(x_dict, demand_dict, p0, pmax, r): new_rate {} for (i, t), var in x_dict.items(): qty var.varValue if qty is None or qty 0: continue n_ops demand_dict[t] / qty # 学习曲线公式 rate_after pmax[(i, t)] - (pmax[(i, t)] - p0[(i, t)]) * np.exp(-r[(i, t)] * n_ops) new_rate[(i, t)] min(rate_after, pmax[(i, t)]) return new_rate rate {(i, t): p0[(i, t)] for (i, t) in price} for iteration in range(15): model, x build_mip(rate) model.solve() if model.status ! 1: print(f第{iteration1}次迭代无解退出) break new_rate simulate_learning(x, demand, p0, pmax, r) change sum(abs(new_rate[k] - rate[k]) for k in new_rate) rate.update(new_rate) print(f迭代{iteration1}: 成本{value(model.objective):.2f}, 成功率变化总量{change:.4f}) if change 1e-3: break print(\n最终订购方案) for (i, t), var in x.items(): if var.varValue and var.varValue 0: print(f供应商{i}, 手术{t}, 订购{int(var.varValue)}台, 学习后成功率{rate[(i,t)]:.4f})代码逻辑说明pulp的LpVariable创建全集lpSum用于建立线性表达式成功率约束中rate[(i,t)]是迭代前的估计每次迭代被新值覆盖simulate_learning假设每台机器人服务demand[t]/qty例手术算出学习后成功率。注意如果某供应商在上一轮没有被采购new_rate里就没有它的项rate.update会保留旧值这符合“没买就不用学”的直觉。参数设置上15次迭代上限和1e-3收敛阈值不是固定值。如果你的数据学习率r很小需要更多迭代反之如果r大于0.1两三步就收敛。min_rate是硬约束如果题目不允许权重平均你要改成“每台机器人的成功率都不得低于阈值”那约束就变成rate[(i,t)] min_rate * x[(i,t)]的分支约束代码相应调整。4.3 怎么读输出成本、购买量、学习后成功率运行上面代码你会看到类似输出迭代1: 成本2850.00, 成功率变化总量0.5160 迭代2: 成本2850.00, 成功率变化总量0.0230 ... 最终订购方案 供应商S2, 手术T1, 订购12台, 学习后成功率0.9634 供应商S2, 手术T3, 订购10台, 学习后成功率0.9477 ...这时候需要注意如果学习后成功率仍然低于min_rate说明数据里根本没有满足条件的组合得调整min_rate或者增加“允许供应商以更高价格提供更高初值成功率”的选项。竞赛里一般不会出现无解但万一出现请优先检查是否把不适配的数据漏过滤了。另一个读输出技巧把“成功率变化总量”打印出来如果连续两次都很小说明学习优化已经稳定。如果这个值震荡不降说明你的迭代逻辑有循环依赖比如某次订购量变小导致每台手术例数变多学习后成功率上升下一轮模型因为成功率变高而减少采购又导致每台手术例数更多……。这种震荡是小规模数据的正常现象缓解办法是把n_ops用上一轮订购量计算但更新rate时乘以一个阻尼系数0.5让成功率更平滑地变化。5. 避坑与常见问题5个数模竞赛里反复出现的建模失误5.1 现象模型解出来成功率约束明明满足但论文里算的正确率对不上原因混合整数规划里用“加权平均成功率”做线性近似但论文后续模拟时用了单台机器人的独立成功率。求得的订购方案保证的是整体加权平均达标不代表每一台机器人都达标。比如按平均算95%可能其中一台是88%另一台是100%平均满足了但实际使用中那台88%的机器会被单独使用某些手术整体成功率就不达标。解决在约束里增加“每台机器人最低成功率线”条件或者用最坏情况成功率做保守约束。如果题目明确“每台机器人都必须达到某值”那么把约束改成for (i, t) in x: if x[(i, t)].varValue: # 只在订购时生效 model rate[(i, t)] min_rate # 注意这在迭代中才有效但这不是线性约束不能直接写进去。一个变通做法是在目标函数里加一个惩罚项对rate min_rate的采购量施加高额惩罚成本。这样模型会在成本和达标之间自行权衡。5.2 现象代码复现时用scipy.optimize.linprog跑整数规划直接报错原因linprog默认只支持连续性变量很多参赛者不知道它不支持整数变量于是把整数约束强制设为bounds(0, None)导致结果出现小数台数成本和成功率全乱套。这是“代码复现”时最常见的翻车点。解决要么换pulp/ortools这类支持整数的库要么用scipy.optimize.milpSciPy 1.9 才有。如果你只能用linprog至少把x的整数性放在论文里说明是“松弛后取整”。但取整可能破坏需求约束我建议直接使用pulp它对竞赛选手极友好安装只需pip install pulp。5.3 现象学习曲线参数敏感稍微改一下初始值最优订购方案就完全不同原因题目往往把学习率r、上限pmax给得含糊甚至只给一个“学习因子”。参赛者凭感觉定参后输出结果随机性大看起来像个黑匣子。这其实是所有“学习优化”模型的通病——参数不是从数据里估计出来的而是拍的。解决做灵敏度分析把r和pmax分别设成低、中、高三档跑九组实验观察订购方案是否稳定。如果某组数据下方案从便宜货跳到了高价货说明模型对参数极度敏感论文里要重点讨论这个现象同时给出“稳健方案”——即九组实验中都出现频率最高的订购组合。更进一步可以从题目给的少量历史数据里用最大似然估计学习参数而不是拍脑袋。5.4 现象迭代10次后成本还没稳定越迭代越贵原因前面说的“成功率上升 → 可减少订购 → 每台手术例数变多 → 成功率继续上升”的反馈环在整数规划下会产生振荡。尤其当需求很小、订购量只有个位数时把8台减到7台会让每台手术例数跳变成功率增量也跳变。解决给rate更新加阻尼同时把手术例数计算改为“四舍五入到底”alpha 0.5 # 阻尼系数 for k in new_rate: rate[k] alpha * new_rate[k] (1 - alpha) * rate.get(k, new_rate[k])这相当于成功率不会因为一次迭代就大幅跳变收敛性会好很多。另一个办法是固定前两轮的订购量只更新一次成功率作为最终方案把两次迭代之间的差写成“学习优化带来的成功率提升”放在论文里反而更有说服力。5.5 现象代码诊断插件报错说LpVariable没有varValue明明求解成功了原因var.varValue只有在model.solve()之后才会被赋值如果求解状态不是OptimalvarValue可能仍是None。你在打印时直接.varValue取值如果模型无解就会拿到None。这是新手最容易踩的坑。解决取值前先判断if model.status 1:pulp的状态码中1代表最优。另外不要在迭代中频繁调用build_mip后忘记重新求解否则上一次的varValue会残留。我习惯在每个循环里重写model和x而不是复用旧对象这样避免变量混淆。若需要调试打印model.status的LpStatus[model.status]看具体状态。6. 验证与进阶用历史数据回测你的订购策略写完求解方案只算完成了70%剩下30%是验证策略在历史数据上是否真的有效。竞赛题一般会给一段历史手术记录格式可能是“日期、供应商、手术类型、成功与否”。你可以用它来模拟“按你的订购方案采购 → 执行手术 → 更新贝叶斯参数 → 下个月重新订购”的滚动过程这比一次买齐所有机器人更贴近实际。回测代码的思路如下# 历史手术记录示意日期、供应商、手术类型、是否成功 import random random.seed(42) history [] for day in range(90): for _ in range(4): surgery random.choice([T1, T2, T3]) supplier random.choice([S1, S2, S3]) success random.random() 0.9 history.append((day, supplier, surgery, 1 if success else 0)) # 按天滚动每30天根据当前Beta分布重算一次订购 from collections import defaultdict beta_params defaultdict(lambda: [1, 1]) results [] for day in range(0, 90, 30): month_data [h for h in history if h[0] day and h[0] day30] # 用本月成功率更新Beta for _, supplier, surgery, ok in month_data: a, b beta_params[(supplier, surgery)] if ok: a 1 else: b 1 beta_params[(supplier, surgery)] [a, b] # 使用当前Beta期望作为rate重新求解订购模型 current_rate {(i, t): a/(ab) for (i, t), (a, b) in beta_params.items()} # 调用build_mip并求解记录本阶段成本和订购量 model, x build_mip(current_rate) model.solve() total_cost value(model.objective) results.append((day, total_cost, current_rate)) print(f第{day}天后: 阶段成本{total_cost:.2f}, 当前成功率{ {k: round(v,3) for k,v in current_rate.items()} })这段回测框架和“python量化交易策略代码”里的回测逻辑很像按周期调仓用过去的数据更新对未来的预期再重新分配资源。它能直观地展示“学习优化”的价值——随着Beta参数越积越多成功率估计越来越准后几个阶段的订购成本通常会低于第一阶段因为不再需要为低成功率“过拟合”地多买冗余机器人。进阶方向有两个一是把“每次采购都有固定订购成本”加进模型让模型自动决定是否该补货二是把学习曲线从Beta-二项分布推广到时间衰减模型比如最近半年的手术记录权重更高因为设备老化会导致成功率下降。一个可用的技巧是给Beta更新时的增量加一个衰减系数alpha decay * alpha success这样老经验会缓慢淡出我更倾向用这个方法处理长期数据。回测完成后我会把各阶段的成本、成功率打印成一张表再和“不做学习优化、只看初始成功率一次买齐”的方案对比。两者成本差就是学习优化的货币化价值。这个对比写进论文比任何理论推导都管用。最后提个习惯我每次跑完这套流程都会把seed固定住确保结果可复现。数据处理脚本、MIP模型、回测脚本分开存不要揉在一个文件里。题目变了改数据表比改代码快得多。希望这些思路和代码能帮你在类似问题上少走几步弯路。本文还有配套的精品资源点击获取
