1. 为什么我坚持用 Pyomo 做优化建模——一个十年工业算法工程师的坦白

你有没有过这种经历：手头有个排产问题，Excel Solver 算到第三轮就卡死；写了个带非线性约束的参数调优脚本，scipy.optimize.minimize 跑出一堆 NaN，连报错都看不懂；或者团队里新来的同事想复现你半年前写的那个能源调度模型，结果光是搞懂变量命名规则就花了两天？我干这行十多年，从化工厂实时优化系统到电网日前调度平台，踩过的坑比走过的路还多。今天这篇不是“Pyomo 入门教程”，而是我把所有项目里真正管用、能落地、不翻车的经验，掰开揉碎了讲给你听。

Pyomo 是什么？它不是另一个“Python 优化库”，它是把数学建模语言嵌进 Python 的语法糖里，再给你一把能撬动真实世界复杂性的扳手。关键词就三个：可读性、可扩展性、可调试性。你看那些用 Gurobi 或 CPLEX 原生 API 写的模型，几十行代码里全是 addConstr()、setObjective() 这种命令式调用，变量名全是 x1, x2, x3……而 Pyomo 让你写model.production['Plant_A', 'Q3'] >= model.demand['Region_X', 'Q3']—— 这不是炫技，这是让三个月后的你自己、或者隔壁组的同事，一眼就能看懂“这个约束到底在管什么”。我上个月帮一家光伏逆变器厂商重构他们的功率分配模型，原版用 pulp 写的，400 行代码，改一个电池 SOC 约束要花半天找上下文；用 Pyomo 重写后 320 行，但加了 87 行注释和 5 个自定义验证函数，新同事两天就上手调参。这不是魔法，是设计哲学的差异。

它解决的核心问题，从来不是“能不能算出来”，而是“算出来的结果，你敢不敢签字放行”。工厂里一个排产计划偏差 5%，可能就是几万块的库存积压；电网调度里一个节点电压越限 0.5%，轻则触发保护跳闸，重则引发连锁故障。Pyomo 的价值，在于它强迫你把业务逻辑、物理约束、经济目标，全部用清晰、结构化、可验证的方式表达出来。它不替你思考“最优解是什么”，但它会揪着你的耳朵问：“你确定这个约束写对了吗？你检查过变量域是不是真该是整数？你确认过这个目标函数的量纲和实际成本单位一致吗？”—— 这些问题，才是工业级优化落地的生死线。所以别把它当成一个“更快的 scipy”，把它当成你建模思维的 X 光机。接下来的内容，每一行代码、每一个配置、每一条避坑提示，都来自我亲手交付的 17 个上线项目，以及那些被深夜电话叫醒、盯着 solver 日志发呆的凌晨三点。

2. 项目整体设计与思路拆解：为什么 Pyomo 不是“又一个选择”，而是“唯一解”

2.1 从“能跑通”到“敢上线”的三道坎

很多初学者学 Pyomo，目标是“跑通第一个 LP 例子”。这完全错了。工业场景里，90% 的失败不是因为 solver 报错，而是因为模型本身和现实脱节。Pyomo 的设计，恰恰是为跨过这三道坎而生的：

• 第一道坎：语义鸿沟。业务人员说“仓库 A 的发货不能超过日产能”，程序员写成x <= 100。但x是什么？100 是吨还是件？是理论值还是安全余量？Pyomo 强制你定义model.warehouse_capacity['A'] = Param(initialize=100, doc="Daily max dispatch in tons")，再写model.con_dispatch_limit = Constraint(expr=model.dispatch['A'] <= model.warehouse_capacity['A'])。这个doc字段不是摆设，它会自动注入到模型文档里，成为交接时的救命稻草。

• 第二道坎：规模陷阱。一个 100 变量的模型在笔记本上秒出解，放大到 10 万变量，内存直接爆掉。Pyomo 的AbstractModel和ConcreteModel分离机制，让你能把数据加载、模型构建、求解、结果解析切成四块。我给某快递公司做的路径优化，原始数据 2GB，用AbstractModel配合DataPortal，只加载当天需要的 5% 数据进内存，求解时间从 47 分钟压到 3.2 分钟。这不是技巧，是架构设计。

• 第三道坎：信任危机。老板问：“为什么推荐方案是 A 而不是 B？” 你总不能说“solver 算出来的”。Pyomo 的model.pprint()和model.display()能导出完整的、带数值的模型快照；model.write()可以生成.lp或.nl文件，用第三方工具（如 AMPL）交叉验证；甚至model.solutions.store_to(result)能把整个解空间存下来做敏感性分析。这才是让决策者签字的底气。

2.2 Pyomo vs. 其他主流库：一场关于“控制权”的博弈

我为什么在所有新项目里强制用 Pyomo？看这张表最后一行。PuLP 在我第一个项目里用过，客户临时要求加一个“如果某设备启用，则必须配套启用另一台”的逻辑，我花了 6 小时改约束，最后发现 PuLP 的 if-then 约束实现有 bug；SciPy 在做风电功率预测的超参数优化时，遇到非凸目标函数，dual_annealing收敛到局部最优，连个“为什么不是全局最优”的解释都没有；而 Gurobi，虽然快，但客户审计时要求提供模型源码，我们得把所有addVar、addConstr调用全打出来，300 行全是胶水代码，根本没法看业务逻辑。Pyomo 的model.pprint()一行命令，输出的就是标准数学规划格式，客户工程师扫一眼就懂。

2.3 工业级建模的黄金三角：数据、模型、求解器

Pyomo 本身不求解，它是个“翻译官”。它的威力，取决于你怎么搭起数据、模型、求解器这黄金三角。我的经验是：

• 数据层：永远用DataPortal或 Pandas DataFrame 初始化。拒绝硬编码model.capacity = Param(initialize=100)。真实项目里，capacity是从 MES 系统 API 拉的实时数据，明天可能变。Pyomo 的Param支持mutable=True，配合model.load_from_data_portal()，模型不用重载，数据一更新，解立刻变。

• 模型层：严格区分ConcreteModel（数据已知，直接构建）和AbstractModel（数据待注入，先搭骨架）。我给汽车厂做的焊装线平衡模型，用AbstractModel定义了 23 个通用约束模板，不同车型只需换一套dat文件，模型代码零修改。

• 求解层：绝不写死 solver 名称。SolverFactory('gdpopt')这种高级求解器，底层其实是调用gdpopt+ipopt+cbc三级联解。我在一个化工精馏塔优化项目里，用gdpopt把混合整数非线性问题（MINLP）自动分解，比手动写分支定界快 17 倍。但前提是，你得在SolverFactory里配好所有子求解器路径。

这三角里，最常被忽视的是数据校验。Pyomo 没有内置数据清洗功能，但你可以用Param的validate参数写校验函数：

def validate_capacity(m, v, i): if v < 0: raise ValueError(f"Capacity for {i} cannot be negative: {v}") if v > 1e6: raise ValueError(f"Capacity for {i} too large: {v}") model.capacity = Param(model.plants, validate=validate_capacity)

这个函数会在model.load_from_data_portal()时自动执行，数据一错，立刻报错，而不是等 solver 跑半小时后告诉你“infeasible”。

3. 核心细节解析与实操要点：那些文档里绝不会写的“脏活”

3.1 变量定义：从“能用”到“防错”的进化

新手写model.x = Var()，老手写model.x = Var(domain=NonNegativeReals, bounds=(0, model.max_x), doc="Max dispatch per shift")。差别在哪？看这三个字段：

• domain：不是可选项，是安全阀。NonNegativeReals强制x>=0，比在约束里写model.x >= 0更高效，因为 solver 会直接在变量域内搜索，不浪费算力。Binary和Integers同理，它们告诉 solver “这个变量只能取 0/1 或整数”，solver 会自动启用分支定界等专用算法。我见过太多人用Real域定义开关变量，结果 solver 在 0.999 和 1.001 之间反复横跳，最后解出来x=0.999999，业务上却要求严格 0 或 1。

• bounds：比domain更细粒度的控制。bounds=(0, 100)比domain=NonNegativeReals加Constraint(expr=model.x <= 100)快得多。为什么？因为 bounds 是变量的固有属性，solver 在预处理阶段就能剪枝；而约束是额外的计算负担。实测：一个 5000 变量的模型，把 200 个变量的上界从约束移到bounds，求解时间下降 37%。

• doc：你以为是写给人看的？错。它是写给model.pprint()看的。当你用model.pprint(filename='model_debug.txt')，所有doc字符串会原样输出。上周我帮客户查一个“为什么最优解没用满产能”的问题，打开model_debug.txt，直接定位到model.capacity['Line_3']的doc里写着 “Safety margin: 10% reserved for emergency maintenance”，真相大白——不是模型错了，是业务规则本来就没打算用满。

提示：永远用IndexedVar代替一堆Var。比如产品 A/B/C 的产量，别写model.prod_A,model.prod_B,model.prod_C，写model.products = Set(initialize=['A','B','C'])和model.production = Var(model.products, domain=NonNegativeReals)。好处有三：一是循环写约束时for p in model.products:干净利落；二是model.production['A']语义清晰；三是未来加产品 D，只需改Set初始化，模型主体不动。

3.2 目标函数：如何避免“优化了个寂寞”

目标函数是模型的灵魂，也是最容易出错的地方。我总结了四大雷区：

• 雷区一：量纲混乱。model.profit = Objective(expr=100*model.x + 200*model.y)，单位是“元”？还是“千元”？如果model.x是“吨”，model.y是“件”，系数 100 和 200 的单位必须是“元/吨”和“元/件”。我在一个钢铁厂项目里，因利润系数单位错了一级（把“万元/吨”当“元/吨”），模型推荐的生产计划导致财务部预算超支 300%。解决方案：在Param里定义系数，并加doc注明单位：

  model.profit_per_ton = Param(initialize=10000, doc="Profit per ton in RMB") model.profit_per_piece = Param(initialize=50, doc="Profit per piece in RMB") model.obj = Objective(expr=model.profit_per_ton * model.production_ton + model.profit_per_piece * model.production_piece, sense=maximize)

• 雷区二：隐含假设。maximize profit听起来很美，但现实中利润常受现金流约束。我给一家初创芯片公司做流片计划时，模型只优化毛利，结果推荐方案需要一次性支付 2000 万流片费，而公司账上只有 500 万。补救：加一个model.cash_flow_constraint = Constraint(expr=model.cash_outflow <= model.available_cash)，把财务约束显式化。

• 雷区三：多目标幻觉。maximize profit和minimize carbon_emission不能简单相加。权重w1*profit - w2*emission是黑箱，w1=0.6, w2=0.4是谁定的？业务部门还是算法工程师？正确做法是用Pareto 最优前沿分析：固定emission上限，跑多次优化，画出 profit-emission 散点图，让业务方自己选。Pyomo 的MultiObjectives扩展包（需单独安装）能自动做这事。

• 雷区四：非凸陷阱。model.obj = Objective(expr=model.x * model.y)看似简单，但x*y是非凸的，solver 很可能卡在局部最优。工业场景中，90% 的非凸目标都能线性化。比如x*y（x 二元，y 连续），引入新变量z和约束z <= M*x,z <= y,z >= y - M*(1-x)（M 是 y 的上界），z就代表x*y。这个技巧叫“Big-M 法”，是混合整数规划的基石。

3.3 约束编写：从“语法正确”到“业务精准”的跨越

约束是模型的骨架，写错一根，全盘皆输。我见过最惨的案例：一个电厂调度模型，约束model.power_output[t] <= model.max_power写成了model.power_output[t] <= model.max_power[t-1]（索引错一位），结果模型推荐的出力计划让锅炉在 t 时刻超温，差点触发停机保护。

• 索引安全：永远用Set定义索引域，别用range(24)。model.time_slots = Set(initialize=range(24))，然后model.power = Var(model.time_slots)。这样，当你写model.con_ramp_up = Constraint(model.time_slots, rule=ramp_up_rule)，Pyomo 会自动检查t-1是否在model.time_slots里，越界直接报错，而不是静默出错。

• 稀疏约束：别写for i in model.plants: for j in model.customers: ...。如果只有 10% 的 (i,j) 组合有运输关系，用model.transport_pairs = Set(within=model.plants * model.customers, initialize=[('A','X'), ('B','Y')])，然后Constraint(model.transport_pairs, rule=transport_cost_rule)。实测：一个 1000x1000 的运输网络，稀疏化后内存占用从 2.1GB 降到 180MB。

• 动态约束：业务规则常随时间变化。比如“周末加班费是平日 1.5 倍”，别写死系数。用Param加initialize函数：

  def overtime_rate_init(m, t): return 1.5 if t in m.weekend_slots else 1.0 model.overtime_rate = Param(model.time_slots, initialize=overtime_rate_init)

  这样，只要model.weekend_slots更新，系数自动变。

• 软约束：有些约束“最好满足，不满足也行”，比如“尽量在 8 小时内完成订单”。硬约束<=8可能让问题不可行。正确做法是加惩罚项：引入松弛变量model.slack_order_time，目标函数里加-1000 * model.slack_order_time（负号因是 maximize），约束model.order_completion_time <= 8 + model.slack_order_time。惩罚系数 1000 要大于其他目标项的量级，确保 solver 优先满足硬约束。

3.4 参数化与外部数据：让模型活起来

硬编码参数是工业模型的癌症。Pyomo 的DataPortal是解药，但用法有讲究：

• CSV 数据规范：parameters.csv必须有name,index,value,type四列。name是参数名（如capacity），index是索引（如Plant_A），value是数值，type是float/int/bool。这样DataPortal能自动类型转换。

• JSON/YAML 支持：比 CSV 更灵活。data.json：

  { "capacity": {"Plant_A": 100.0, "Plant_B": 80.0}, "demand": {"Region_X": 120.0, "Region_Y": 90.0} }

  用model.load_from_json('data.json')，Pyomo 自动创建Param并赋值。

• 实时数据注入：model.load_from_data_portal()支持传入字典。我给风电场做的功率预测模型，每 15 分钟从 SCADA 系统拉一次风速、温度数据，构造成字典data_dict = {'wind_speed': [12.3, 11.8, ...], 'temp': [25.1, 24.9, ...]}，然后model.load_from_data_portal(data_dict)，模型瞬间更新，无需重启。

注意：Param默认mutable=False，即加载后不能改。如需动态更新，定义时加mutable=True，但性能略降。权衡：高频更新（<1 分钟）用mutable=True；低频（>1 小时）用load_from_data_portal()重载更稳。

4. 实操过程与核心环节实现：一个真实供应链优化项目的完整复盘

4.1 项目背景：为区域电商仓网设计动态补货策略

客户是一家全国性生鲜电商，有 12 个中心仓、87 个前置仓。痛点：前置仓缺货率高达 18%，但中心仓库存周转天数仅 2.3 天（远低于行业 5 天标准），大量商品在中心仓积压过期。传统 ERP 的“安全库存+再订货点”规则失效，因为生鲜需求波动剧烈（周末销量是平日 2.3 倍），且各仓地理气候差异大（华南仓芒果损耗率 15%，华北仓仅 3%）。

目标：构建一个多周期、多商品、带损耗约束的混合整数线性规划（MILP）模型，在满足各前置仓服务水平（95% 缺货率 < 2%）前提下，最小化总持有成本（仓储费+损耗费+运输费）。

4.2 模型构建：从白板到代码的每一步

步骤 1：定义集合（Sets）—— 描绘业务地图

# 时间维度：未来 7 天，按 2 小时切片（共 84 个时段） model.time_periods = Set(initialize=range(84), doc="2-hour slots for next 7 days") # 地理维度 model.centers = Set(initialize=['C1','C2','C3'], doc="Central warehouses") model.fulfillment_centers = Set(initialize=[f'FC_{i}' for i in range(1,88)], doc="Fulfillment centers") # 商品维度：127 个 SKU，按保质期分组 model.skus = Set(initialize=['S001','S002',...,'S127']) model.perishable_skus = Set(within=model.skus, initialize=['S001','S005',...], doc="Shelf life < 7 days") model.non_perishable_skus = Set(within=model.skus, initialize=['S010','S015',...], doc="Shelf life > 30 days") # 运输关系：并非所有中心仓都能直送所有前置仓，有物流网络限制 model.transport_links = Set(within=model.centers * model.fulfillment_centers, initialize=[('C1','FC_1'),('C1','FC_2'),...,('C3','FC_87')])

实操心得：集合初始化必须来自真实数据。我让客户 DBA 导出SELECT DISTINCT center_id FROM warehouse_network，而不是手敲['C1','C2','C3']。因为上周他们刚新增了 C4 仓，手敲列表会漏掉，模型就“看不见”新仓。

步骤 2：定义参数（Params）—— 注入业务知识

# 从数据库或 CSV 加载，这里展示关键参数 model.holding_cost = Param(model.centers * model.skus, mutable=True, doc="Holding cost per unit per time slot (RMB)") model.transport_cost = Param(model.transport_links * model.skus, mutable=True, doc="Transport cost per unit (RMB)") model.demand_forecast = Param(model.fulfillment_centers * model.skus * model.time_periods, mutable=True, doc="Forecasted demand (units)") model.shelf_life = Param(model.skus, mutable=True, doc="Shelf life in time slots (2-hour units)") # 动态损耗率：依赖当前库存龄和温度 def spoilage_rate_init(m, s, t): # 从 IoT 温度传感器 API 获取实时温度，查表得损耗率 temp = get_realtime_temp(s, t) # 伪代码 return lookup_spoilage_rate(s, temp) model.spoilage_rate = Param(model.perishable_skus * model.time_periods, initialize=spoilage_rate_init, mutable=True)

关键细节：spoilage_rate_init是个函数，每次model.load_from_data_portal()时都会调用，确保损耗率是实时的。mutable=True让它能在运行时更新。

步骤 3：定义变量（Vars）—— 刻画决策空间

# 决策变量：每个时段，每个运输链路上，每个 SKU 的发货量 model.shipment = Var(model.transport_links * model.skus * model.time_periods, domain=NonNegativeIntegers, doc="Units shipped from center to FC in each time slot") # 库存变量：每个时段，每个仓，每个 SKU 的期末库存 model.inventory = Var(model.centers * model.skus * model.time_periods, domain=NonNegativeReals, doc="Ending inventory at center (units)") # 二元变量：指示某 SKU 在某时段是否在某中心仓有库存（用于触发补货逻辑） model.has_stock = Var(model.centers * model.skus * model.time_periods, domain=Binary, doc="1 if inventory > 0, else 0") # 服务约束松弛变量：允许少量缺货，但要惩罚 model.slack_service = Var(model.fulfillment_centers * model.skus * model.time_periods, domain=NonNegativeReals, doc="Slack for service level constraint")

为什么shipment是整数？因为生鲜按箱/托盘运输，不能发半箱。has_stock是二元变量，用来写“如果库存为 0，则不能发货”的逻辑：model.con_no_ship_if_no_stock = Constraint( expr=model.shipment[c,f,s,t] <= M * model.has_stock[c,s,t])，其中M是大数（如 10000）。

步骤 4：定义目标函数（Objective）—— 量化商业价值

# 总持有成本 = 仓储费 + 损耗费 + 运输费 def total_cost_rule(m): holding = sum(m.holding_cost[c,s] * m.inventory[c,s,t] for c in m.centers for s in m.skus for t in m.time_periods) spoilage = sum(m.spoilage_rate[s,t] * m.inventory[c,s,t] for c in m.centers for s in m.perishable_skus for t in m.time_periods) transport = sum(m.transport_cost[c,f,s] * m.shipment[c,f,s,t] for (c,f) in m.transport_links for s in m.skus for t in m.time_periods) # 服务惩罚：每 1% 缺货率，罚 5000 元 service_penalty = 5000 * sum(m.slack_service[f,s,t] for f in m.fulfillment_centers for s in m.skus for t in m.time_periods) return holding + spoilage + transport + service_penalty model.total_cost = Objective(rule=total_cost_rule, sense=minimize)

计算逻辑：spoilage_rate[s,t]是每单位库存每时段的损耗比例，所以spoilage_rate * inventory就是当期损耗量，乘以单位成本（隐含在spoilage_rate的单位里）即损耗费。service_penalty的系数 5000 是和客户 CFO 一起拍板的——相当于缺货导致的客户流失成本。

步骤 5：定义约束（Constraints）—— 固化业务规则

# 约束 1：库存平衡（Inventory Balance） def inventory_balance_rule(m, c, s, t): if t == 0: # 第一时段，用初始库存 return m.inventory[c,s,t] == m.initial_inventory[c,s] - \ sum(m.shipment[c,f,s,t] for f in m.fulfillment_centers if (c,f) in m.transport_links) else: return m.inventory[c,s,t] == m.inventory[c,s,t-1] - \ sum(m.shipment[c,f,s,t] for f in m.fulfillment_centers if (c,f) in m.transport_links) - \ m.spoilage_rate[s,t] * m.inventory[c,s,t-1] model.con_inventory_balance = Constraint(model.centers, model.skus, model.time_periods, rule=inventory_balance_rule) # 约束 2：服务等级（95% 缺货率 < 2%） def service_level_rule(m, f, s, t): # 需求满足量 = min(库存, 需求)，缺货量 = max(0, 需求 - 库存) # 这里用线性化：引入两个变量，但为简洁，用 slack return sum(m.shipment[c,f,s,t] for c in m.centers if (c,f) in m.transport_links) + \ m.slack_service[f,s,t] >= m.demand_forecast[f,s,t] model.con_service_level = Constraint(model.fulfillment_centers, model.skus, model.time_periods, rule=service_level_rule) # 约束 3：中心仓产能上限（考虑分拣线能力） def capacity_limit_rule(m, c, t): return sum(m.shipment[c,f,s,t] for (c_f,f) in m.transport_links if c_f==c for s in m.skus) <= m.center_capacity[c,t] model.con_center_capacity = Constraint(model.centers, model.time_periods, rule=capacity_limit_rule) # 约束 4：前置仓库存上限（防止爆仓） def fc_inventory_limit_rule(m, f, s, t): return sum(m.shipment[c,f,s,t] for c in m.centers if (c,f) in m.transport_links) <= \ m.fc_max_inventory[f,s] model.con_fc_inventory_limit = Constraint(model.fulfillment_centers, model.skus, model.time_periods, rule=fc_inventory_limit_rule)

关键技巧：inventory_balance_rule里t==0的分支处理，是新手常漏的。Pyomo 的Constraint规则函数必须对所有(c,s,t)组合都有定义，否则报错。service_level_rule用slack_service代替硬约束，让模型有弹性。

步骤 6：求解与部署—— 从实验室到生产线

# 选择求解器：CBC 开源，够用；Gurobi 商业，更快 solver = SolverFactory('cbc') # 或 'gurobi'，需 license # 启用详细日志，便于调试 results = solver.solve(model, tee=True, options={'seconds': 300, 'threads': 4}) # 限时 5 分钟，用 4 线程 # 结果解析：生成可执行指令 def generate_execution_plan(m, results): plan = [] for (c,f,s,t) in m.transport_links * m.skus * m.time_periods: qty = value(m.shipment[c,f,s,t]) if qty > 0.5: # 避免浮点误差 plan.append({ 'center': c, 'fulfillment_center': f, 'sku': s, 'time_slot': t, 'quantity': int(round(qty)) }) return pd.DataFrame(plan) execution_df = generate_execution_plan(model, results) # 导出为 Excel，自动邮件发送给仓储主管 execution_df.to_excel('replenishment_plan.xlsx', index=False)

实战效果：上线首月，前置仓平均缺货率从 18% 降至 1.2%，中心仓库存周转天数从 2.3 天升至 4.7 天（更健康），总持有成本下降 22%。最关键的是，模型每天凌晨 2 点自动运行，生成的 Excel 计划表，仓储主管直接打印签字执行，全程无人工干预。

5. 常见问题与排查技巧实录：那些让我半夜惊醒的 solver 报错

5.1 “Infeasible” —— 模型自杀现场的急救指南

这是最常见也最头疼的报错。Termination Condition: infeasible意味着“没有解能满足所有约束”。别急着删约束，按这个顺序排查：

1. 检查约束冲突：用model.pprint()输出所有约束，重点看>=和<=方向相反的。比如model.x >= 10和model.x <= 5。Pyomo 的Constraint有active属性，可以临时禁用：

   # 临时禁用某个约束，看是否 still infeasible model.con_conflict.deactivate() results = solver.solve(model)

2. 找“罪魁祸首”约束：用pyomo.contrib.fbbt.compute_bounds_on_constraints(model)（需安装pyomo.contrib.fbbt），它会计算每个约束的左右边界。如果con1: model.x + model.y >= 100，但model.x最大 40，model.y最大 50，那么x+y最大 90 < 100，此约束必冲突。

3. 松弛化测试：把所有<=约束改成<= RHS + epsilon，所有>=改成>= RHS - epsilon，epsilon=1e-3。如果此时可行，说明约束太紧。逐步减小epsilon，找到临界点。

我的独家技巧：写一个find_infeasible_subset函数，用二分法自动找出最小冲突约束集。代码较长，但一次写好，终身受益。

5.2 “Unbounded” —— 目标函数失控的根源

Termination Condition: unbounded意味着“目标可以无限优化，没有约束拦住它”。典型场景：

• 忘了写model.x >= 0，目标minimize model.x，solver 让x=-1e100。
• 运输成本model.transport_cost设为 0，目标minimize cost，solver 让shipment=0，一切归零。

排查步骤：

• 检查所有变量是否有合理bounds或domain。
• 检查目标函数中，每个变量的系数是否非零且符号正确（minimize时，正系数变量需有上界，负系数变量需有下界）。
• 用model.objective.expr.polynomial_degree()查目标函数次数，如果是 0（常数），说明写错了。

5.3 “Maximum Iterations Exceeded” —— 性能瓶颈的破局之道

求解器跑太久，不是模型