簡介面向具備線性規(guī)劃基礎(chǔ)與Python編程能力的研究人員與工程師PDF內(nèi)容系統(tǒng)講解如何利用PuLP庫對多目標線性規(guī)劃問題進行建模與求解。內(nèi)容以復(fù)雜約束下的資源分配優(yōu)化為背景完整覆蓋問題初始化、54個決策變量定義、三組目標函數(shù)如最大化RE、最小化Q設(shè)置以及資源、比例、需求滿足、總成本、混合、最大值、變量分解與整數(shù)變量等多類約束的添加并給出求解與結(jié)果輸出的核心步驟。針對多目標處理資料對比分層優(yōu)化法與加權(quán)法提供分層優(yōu)化的可運行代碼及逐步解釋便于讀者遷移到相似優(yōu)化任務(wù)同時討論模型復(fù)雜度提出引入商業(yè)求解器、采用ε-約束法等改進方向。資源包為單個PDF文檔424KB已有84人學習適合用于論文復(fù)現(xiàn)、課程設(shè)計或工程中資源分配方案的驗證。對于需要系統(tǒng)掌握PuLP建模范式的讀者尤為適用。1. 這份資源能解決什么問題54個變量、38條約束的PuLP建模原型做資源分配優(yōu)化的讀者多數(shù)是沖著 PuLP 線性規(guī)劃建模來的但這套代碼最大的價值不是“能跑”而是把論文里的多目標優(yōu)化問題完整落成了可復(fù)現(xiàn)的 Python 模型。它包含 54 個決策變量6 組 × 3 子組 × 3 子子組、三個方向的目標函數(shù)和六類約束幾乎覆蓋了實際工程項目里會遇到的約束類型資源上限、比例限制、需求下限、總成本、混合不等式以及需要線性化才能處理的 max 約束。對正在寫論文復(fù)現(xiàn)、做課程設(shè)計、或者要給企業(yè)資源分配問題搭第一版原型的工程師來說這份代碼是一個可以直接改參數(shù)的起點而不是需要從零推導(dǎo)的數(shù)學題。我拆這份資源時最大的感觸是PuLP 建模真正的門檻不在庫的語法而在約束條件怎么組織、多目標怎么分層、max 這類非線性約束怎么用輔助變量繞過。2. 先把模型拆開54個決策變量與三層目標函數(shù)的表達方式2.1 決策變量怎么命名最不容易抄錯這份資源里的決策變量是 x_ijk含義是“第 i 組、第 j 子組、第 k 子子組的分配量”。如果你的論文里也出現(xiàn)這種三維下標最怕的就是手寫 54 行 LpVariable 聲明抄到后面連索引都亂了。常見做法是先用循環(huán)字典生成變量再用元組當鍵from pulp import * # 創(chuàng)建問題實例先明確主目標是最大化 prob LpProblem(Multi_Objective_Optimization, LpMaximize) # 用字典保存 54 個決策變量 x {} for i in range(1, 7): # 組別 i1..6 for j in range(1, 4): # 子組 j1..3 for k in range(1, 4): # 子子組 k1..3 x[(i, j, k)] LpVariable(fx_{i}{j}{k}, lowBound0)這段代碼的邏輯是用一個三元組 (i, j, k) 作為字典鍵變量名直接拼接成 x_111 這種格式方便求解后輸出時對應(yīng)到論文表格。lowBound0 表示所有決策變量非負這是線性規(guī)劃里最常見的默認假設(shè)。需要特別留意的是假如原問題允許負值或者某些變量必須是整數(shù)就得在這里給每個變量單獨加 catLpInteger 或設(shè)置 lowBound 和 upBound。我一般會先在代碼里跑一個len(x)檢查確認變量數(shù)是 54 而不是 53 或 55這個動作只要花兩秒?yún)s能避免后面所有約束全部錯位。2.2 三個目標的表達RE權(quán)重、Q和輔助變量f/a原問題里有三個方向的目標最大化 RE帶權(quán)收益最小化 Q第 6 組的分配總量以及最小化 fun1f 和 fun2a。這里要看清一個關(guān)鍵點RE 的系數(shù)只跟組別 i 有關(guān)跟 j、k 無關(guān)所以 RE 應(yīng)該是對每個 i 的 9 個變量先求和再乘對應(yīng)權(quán)重。原始代碼里逐項展開寫容易在續(xù)行時丟系數(shù)我在拆解時更推薦先構(gòu)建權(quán)重字典再循環(huán)求和# RE 系數(shù)只與組別 i 相關(guān) re_coeff {1: 0.23, 2: 0.18, 3: 0.20, 4: 0.78, 5: 0.62, 6: 0.00} # 計算 RE對每個組別 i累加該組下 9 個變量的總和再乘權(quán)重 RE lpSum(re_coeff[i] * lpSum(x[(i, j, k)] for j in range(1, 4) for k in range(1, 4)) for i in range(1, 7)) # Q第 6 組的變量總和作為次要目標 Q lpSum(x[(6, j, k)] for j in range(1, 4) for k in range(1, 4)) # 輔助變量 f 和 a用于復(fù)雜約束16不直接進目標 f LpVariable(f, lowBound0) a LpVariable(a, lowBound0) # 第一階段主目標最大化 RE prob.setObjective(RE)邏輯說明用字典存系數(shù)、用 lpSum 做累加比原始寫法少了一堆0.23 * (sum(...))的括號嵌套也不容易在復(fù)制粘貼時漏項。參數(shù)上權(quán)重 {1:0.23, 2:0.18, 3:0.20, 4:0.78, 5:0.62, 6:0.00} 說明第 4 組和第 5 組對 RE 的貢獻度最高第 6 組權(quán)重為 0這意味著最大化 RE 時第 6 組變量在目標函數(shù)里沒有直接激勵它的存在主要是為了滿足約束里的需求下限。f 和 a 在目標函數(shù)里并沒有真正參與它們是給復(fù)雜約束 16 用的輔助變量后面求解時會通過約束把 f 和相關(guān)表達式綁定。2.3 多目標怎么塞進單目標LP分層法與加權(quán)法PuLP 本身不支持多目標函數(shù)這是很多新手卡住的地方。資源里給的策略有兩種我在實際項目里也都是這么做的。分層法的思路是先把主目標 RE 最大化求出最優(yōu)值記為 re_opt再把這個最優(yōu)值當作約束條件放回模型允許一定比例的退化然后切換目標函數(shù)去優(yōu)化次要目標 Q。注意這里第二層求解時要讓 Q 最小化而不能繼續(xù)沿用原來的 LpMaximize 設(shè)置否則跑出來的是最大值——這個細節(jié)我在避坑章會展開。# 第一階段求 RE 最優(yōu)值 prob.solve() re_opt value(prob.objective) # 第二階段允許 RE 有 1% 的下降把 RE 變成約束 prob RE 0.99 * re_opt, RE_epsilon_constraint # 新建一個最小化問題單獨求 Q prob2 LpProblem(Second_Stage, LpMinimize) # 把變量字典復(fù)制進新問題注意PuLP 的變量可以跨問題復(fù)用 prob2.addVariables(x.values()) prob2 Q, Minimize_Q # 把 RE 的 epsilon 約束和其他所有約束重新加進 prob2加權(quán)法則是把多個目標按權(quán)重合成一個prob RE - w * Q但這里的坑是量綱。RE 和 Q 數(shù)值可能差幾十倍w 取 0.1 還是 10 對結(jié)果影響非常大必須結(jié)合具體數(shù)據(jù)量級去調(diào)。分層法的優(yōu)勢在于不需要人為定權(quán)重第一層解就是 RE 的全局最優(yōu)RE 的損減比例0.99也更好理解。我在實際項目里兩種都會跑一遍對比結(jié)果穩(wěn)定性。3. 約束條件分批落地從資源約束到混合約束的可執(zhí)行寫法3.1 資源約束1-6用循環(huán)替代手寫避免系數(shù)抄岔資源約束的部分原始代碼手寫了六個約束每個都是“一個維修周期內(nèi)某種資源的消耗量不能超過上限”。第一個約束是7*(x111x121x131) 2*(x112x122x132) 3*(x113x123x133) 8316。注意到每一行只涉及同一個組別 i 的三個 j、k 組合但系數(shù)跟 j 沒有關(guān)系只跟 k 有關(guān)。這其實可以整理成系數(shù)表用循環(huán)生成比手寫 6 行大括號表達式要清晰得多而且后面如果要改資源上限只需要改一個字典。# 資源系數(shù)resource_coeff[i][k] resource_coeff { 1: {1: 7, 2: 2, 3: 3}, 2: {1: 5, 2: 1.5, 3: 3}, 3: {1: 6, 2: 2, 3: 4}, 4: {1: 5, 2: 2, 3: 3}, 5: {1: 6, 2: 2, 3: 4}, 6: {1: 8, 2: 1.5, 3: 3}, } resource_limit {1: 8316, 2: 13860, 3: 11088, 4: 23100, 5: 27720, 6: 4620} for i in range(1, 7): prob lpSum( resource_coeff[i][k] * x[(i, j, k)] for j in range(1, 4) for k in range(1, 4) ) resource_limit[i], fResource_Constraint_{i}這段代碼把 6 個約束壓縮成 6 行循環(huán)體修改資源上限時只需要改 resource_limit 字典。邏輯上約束名 Resource_Constraint_1 到 6 會被 PuLP 自動編號如果你在求解后要查某個約束的對偶值直接按名字訪問就行。參數(shù)說明第一組資源上限 8316 最小第六組 4620 也小但第四、第五組給了 23100 和 27720說明后兩組是資源投放的重心這也和 RE 系數(shù)權(quán)重 0.78、0.62 是一致的。3.2 比例約束7-15與26-34同一套模板兩套系數(shù)比例約束這塊是資源里最容易讓人心態(tài)崩的地方因為它有九組約束 7-15每組長得很像但系數(shù)和右側(cè)比例完全不同。我拆解時逐條對比了原文代碼發(fā)現(xiàn)這些約束本質(zhì)上是“各組污染率或損耗率的加權(quán)平均不能超過某個上限”右側(cè)乘的是對應(yīng)子組的變量總和。原始代碼里的約束 7 是0.08x111 0.04x211 0.05x311 0.02x411 0.01x511 0.08x611 0.05 * (x111x211...x611)。這里有個值得注意的結(jié)構(gòu)這個約束的左側(cè)和右側(cè)同時包含變量不能簡單地把右側(cè)移到左邊再合并同類項因為左移后 x111 的系數(shù)變成 0.08 - 0.05 0.03實際是把“平均占比不能超過 5%”這種非線性比率約束通過乘以總和一個線性化技巧變成了線性不等式。同類的約束 26-34 也遵循一樣的模式只是系數(shù)不同。寫成代碼可以建立一個 ratio_coeff 字典來復(fù)用模板真正要改的只是每個約束的系數(shù)列表和右側(cè)上限值。# 約束7-15的系數(shù)和上限只列前3組其余按論文表補全 ratio_cfg [ # (右側(cè)比例上限, 系數(shù)元組, 變量位置元組) (0.05, (0.08, 0.04, 0.05, 0.02, 0.01, 0.08), (1,1,1)), (0.04, (0.08, 0.04, 0.05, 0.02, 0.01, 0.08), (2,1,1)), (0.04, (0.08, 0.04, 0.05, 0.02, 0.01, 0.08), (3,1,1)), # ... 約束10-15類似 ] for idx, (limit, coeffs, pos) in enumerate(ratio_cfg, start1): total lpSum(x[(i, *pos)] for i in range(1, 7)) lhs lpSum(coeffs[i-1] * x[(i, *pos)] for i in range(1, 7)) prob lhs limit * total, fRatio_Constraint_{idx 6}這里的循環(huán)變量 pos 是一個二元組代表固定的 (j, k) 位置total 是該位置下 6 個組別的變量總和。為什么要這樣組織因為約束 7-15 每條都只作用于一個 (j, k) 組合把位置元組抽出來循環(huán)生成能避免把 26-34 的系數(shù)抄到 7-15 里去。我在第一次復(fù)現(xiàn)時就把兩組約束的順序搞混了導(dǎo)致結(jié)果無界后來改成這種帶編號的配置結(jié)構(gòu)一眼就能看出每個約束的用途。3.3 需求約束17-25與損耗率(1-loss)的含義需求約束是另一個高頻翻車點因為系數(shù)里出現(xiàn)了 (1-0.08) 這種寫法很多人不理解為什么變量要乘以一個小于 1 的數(shù)再和需求值比較。實際上這表示“實際可用量必須覆蓋需求”也就是每組變量 x 在扣除 8% 或 5% 的損耗后剩余的凈量要滿足對應(yīng)子組的需求下限。原始代碼里第一組需求約束是(1-0.08)x111 (1-0.04)x211 ... 1039這里的 0.08、0.04 就是各組在該子子組下的損耗率1039 是凈需求。代碼里可以定義損耗率矩陣再循環(huán)生成 9 條約束loss_rate { # loss_rate[(i, k)] 表示第 i 組在 k 子子組的損耗率 (1, 1): 0.08, (2, 1): 0.04, (3, 1): 0.05, (4, 1): 0.02, (5, 1): 0.01, (6, 1): 0.08, (1, 2): 0.05, (2, 2): 0.05, (3, 2): 0.04, (4, 2): 0.03, (5, 2): 0.03, (6, 2): 0.05, (1, 3): 0.06, (2, 3): 0.02, (3, 3): 0.04, (4, 3): 0.03, (5, 3): 0.02, (6, 3): 0.06, } demand { (1, 1): 1039, (2, 1): 1732, (3, 1): 1558, (1, 2): 3506, (2, 2): 5844, (3, 2): 4675, (1, 3): 3312, (2, 3): 5979, (3, 3): 4140, } for j in range(1, 4): for k in range(1, 4): prob lpSum((1 - loss_rate[(i, k)]) * x[(i, j, k)] for i in range(1, 7)) \ demand[(j, k)], fDemand_Constraint_{17 (j-1)*3 (k-1)}這個循環(huán)的編號邏輯是按原論文約束 17-25 的順序排的方便和論文表格一一對應(yīng)。結(jié)合常識想一下需求約束是社會需求資源約束是能力上限比例約束是質(zhì)量限制這三類同時存在時模型找的解必須是產(chǎn)能、質(zhì)量、需求三者之間的平衡點。如果跑出來的結(jié)果提示不可行第一步應(yīng)該去檢查是不是需求下限定得比資源上限還高。3.4 總成本約束35與混合約束36-39等式與不等式混排時的檢查順序總成本約束是把 6 個組別不同子子組的單位成本分別乘上對應(yīng)變量累加后不能超過 18300000。這一條在原始代碼里是一行巨大的表達式我拆的時候重新整理成了單位成本字典。混合約束 36-39 則是把 0.5、0.015、0.1 之類的混合系數(shù)組合起來再減掉一個基準值和另一個常數(shù)比較。約束 36 的寫法是0.5*(x111...x611) 0.015*(x112...x612) 0.1*(x113...x613) - 1400 956左邊減 1400 再比較等價于把固定消耗也納入預(yù)算。這組約束在代碼實現(xiàn)上不復(fù)雜但要特別注意的是它們前面有的是 有的是 混在一起時如果統(tǒng)一按 寫模型就直接錯了。我一般把所有約束先按類別注釋分組求解后打印所有約束的 slack檢查哪些約束是緊的slack0哪些約束是松的能快速定位約束方向是不是寫反了。# 單位成本只與組別i和子子組k相關(guān) unit_cost { (1, 1): 1650, (1, 2): 150, (1, 3): 135, (2, 1): 2520, (2, 2): 118, (2, 3): 186, (3, 1): 1560, (3, 2): 130, (3, 3): 80, (4, 1): 2315, (4, 2): 188, (4, 3): 205, (5, 1): 2200, (5, 2): 128, (5, 3): 132, (6, 1): 2006, (6, 2): 110, (6, 3): 62, } prob lpSum(unit_cost[(i, k)] * x[(i, j, k)] for i in range(1, 7) for j in range(1, 4) for k in range(1, 4)) 18300000, \ Total_Cost_Constraint_35混排檢查有一個很實用的習慣用一個列表把約束名、方向、右側(cè)常數(shù)集中管理跑完一遍后統(tǒng)一切換方向做 sanity check。比如成本約束從 改成 模型目標值一定會漲上去如果沒漲說明約束要么根本沒生效要么變量單位不對這個現(xiàn)象比任何調(diào)試工具都管用。4. 避坑手冊PuLP多目標求解的五個翻車點與排查方法4.1 第二層目標還在LpMaximize里求Q根本不會變小現(xiàn)象第一層求完 RE 之后把目標函數(shù)切到 Q再調(diào) prob.solve()輸出的 Q 值居然比第一層的還大完全不是期望的最小化結(jié)果。 原因LpProblem 在初始化時被設(shè)置成了 LpMaximize整個問題的“sense”就固定了。當你 prob.setObjective(Q) 再 solve()PuLP 仍然按最大化方向求解求的是 Q 的最大值而不是最小值。 解決第二層新建一個 LpMinimize 的問題實例把變量、約束全部加進去再求解。我復(fù)現(xiàn)時的做法是寫一個 build_problem(prob, sense) 的函數(shù)把約束函數(shù)抽出來共用兩階段各調(diào)一次避免手寫兩套重復(fù)代碼。4.2 手寫RE表達式括號錯位0.62系數(shù)被截斷之后發(fā)生什么現(xiàn)象原始代碼里 RE 的表達式鏈條特別長復(fù)制粘貼時容易出現(xiàn)括號對不齊續(xù)行末尾反斜杠后多了一個空格或者某一行系數(shù) 0.62 丟了。跑起來不報錯但目標函數(shù)數(shù)值比論文結(jié)果低一大截。 原因PuLP 對表達式是按語法解析的括號沒閉合或續(xù)行符丟失時后面接的內(nèi)容可能被當成新表達式或者直接變成 0 系數(shù)導(dǎo)致 RE 項被靜默丟棄。線性規(guī)劃求解器不會提醒你目標函數(shù)里有變量沒參與。 解決改用字典系數(shù) lpSum 循環(huán)的結(jié)構(gòu)把權(quán)重、索引、值分離跑完第一層后手動打印value(prob.objective)和論文基準值對比。我給自己定的規(guī)矩是任何超過 10 項的目標函數(shù)都不手寫一律先存成 dict再lpSum(coeff[v] * v for v ...)。這一步能擋掉九成以上的表達式 bug。4.3 max約束輔助變量不要把它寫進目標函數(shù)現(xiàn)象加了 max 約束線性化的輔助變量 max_ijk 之后模型解出來的變量分布變得很奇怪某些組別被無意義壓得很低RE 也比預(yù)期低不少。 原因max 約束max{...} 13線性化需要引入一個新變量 z并添加 z 每一項、z 13 這幾條約束。這里的 z 只是傳遞“最大值不能超過 13”這個信息它不該進入目標函數(shù)。但如果寫代碼時圖省事把 z 順手加進了prob RE - 0.1 * sum(max_vars)求解器就有動力把 z 壓到比實際最大值更低從而多余地壓縮了所有相關(guān)變量。 解決輔助變量只出現(xiàn)在約束里不參與任何目標函數(shù)表達式。如果你需要確認哪條 max 約束是緊的求解后打印輔助變量值和對應(yīng)各項的值對比松緊一眼就能看出來。這個坑很有迷惑性因為模型仍然有解只是解被“人為扭曲”了。4.4 無界或不可行先查LpStatus再查對偶值現(xiàn)象模型跑完沒有報錯但 print 出來的變量全為 0或者 LpStatus 顯示 Infeasible然后一堆人開始懷疑是不是 PuLP 裝錯了。 原因不可行多數(shù)時候是約束方向?qū)懛椿蛴覀?cè)常數(shù)抄錯比如把 8326抄成了 8326或者需求下限 1039 抄成了 10390無界則通常是目標函數(shù)里有變量但缺少上界約束。PuLP 默認的 CBC 求解器在遇到不可行時并不會直接告訴你哪條約束出問題。 解決第一件事是檢查prob.status如果是 1 繼續(xù)往下走是 0 或者 -1 就要逐類排查。我一般會寫一小段代碼遍歷prob.constraints.items()打印每個約束的 slack 和雙變量對偶值哪個約束的 slack 是接近 0 且對偶值很大哪個就是模型里最“緊”的環(huán)節(jié)。如果還是定位不到就把約束按類別注釋掉一部分重跑用二分法找是哪一類約束引起的矛盾。這個方法笨但絕對有效。4.5 數(shù)值噪聲當非零解輸出前先做epsilon過濾現(xiàn)象打印“所有非零變量”時出現(xiàn)了 x_111 1e-9、x_222 2e-8 這種近乎為 0 的值看起來好像算出了一堆變量實際模型解幾乎是空的。 原因浮點求解器在迭代收斂時變量值會有數(shù)值噪聲CBC 尤其常見。直接把if v.varValue 0作為篩選條件會把 1e-9 這種噪聲當成有意義的分配量影響后續(xù)統(tǒng)計 Q 和總成本時出現(xiàn)幾行毫無意義的“微量分配”。 解決判斷閾值設(shè)到 1e-5 或 1e-6輸出時統(tǒng)一 round 到四位小數(shù)同時要留意varValue本身可能為 None變量沒參與任何約束時直接比較會報錯。我把篩選條件寫成if v.varValue is not None and v.varValue 1e-5這就不會再被數(shù)值噪聲騙了。5. 驗證與進階求解完成后怎么確認解可靠、怎么掃出帕累托前沿5.1 四個檢查動作狀態(tài)、目標值、冗余約束、對偶檢驗求解完不要直接拿變量去寫論文我習慣按順序做四個檢查。第一打印LpStatus[prob.status]必須是 Optimal 才繼續(xù)否則返回去改約束。第二核對主目標 RE 的值是否落在論文給出的量級范圍比如論文基準 RE 在某個量級你算出來差 10 倍那一定是權(quán)重或系數(shù)單位錯了。第三打印所有約束的 slack找出冗余約束——有些約束在最優(yōu)解下 slack 遠大于 0說明它根本沒起約束作用這類約束可以在后續(xù)做靈敏度分析時排除減少模型復(fù)雜度。第四檢查對偶值對偶值大說明對應(yīng)約束是資源瓶頸比如 Resource_Constraint_1 的對偶值很高意味著增加 8316 那個上限能顯著提升 RE這是給企業(yè)的資源投放建議里最值錢的信息。# 常見做法求解后統(tǒng)一輸出狀態(tài)、目標值和關(guān)鍵約束的松弛量 print(Status:, LpStatus[prob.status]) print(RE , value(prob.objective)) for name, constraint in prob.constraints.items(): slack constraint.slack if slack is not None and abs(slack) 1e-6: print(fTight constraint: {name}, dual {constraint.pi})這段代碼里constraint.pi是對偶值它只有在約束是緊的時候才有意義。如果一個約束 slack 特別大pi 會是 0說明當前資源上限再怎么加也不會提升目標值。做項目匯報時把這三個數(shù)亮出來比貼一整頁求解日志更有說服力。5.2 用ε約束掃描RE閾值得到多組備選方案分層優(yōu)化的第二階段RE 的退化比例是 0.99也就是允許 RE 下降 1% 去換 Q 的降低。但實際決策場景里你可能想知道“RE 每降低 1%Q 能降多少”這就需要對 ε 做參數(shù)掃描。做法是循環(huán)調(diào)整 RE 的保留比例從 0.99 逐步降到 0.95每跑一次記錄一組 (RE, Q)得到一系列備選方案。把這幾組數(shù)據(jù)畫出來就是近似的帕累托前沿決策者可以直觀地看到收益和成本之間的權(quán)衡關(guān)系。從那次以后我每次用這個模型都會強制跑一遍閾值掃描哪怕最后只取其中一組解也要讓決策者看到還有哪些備選方案。參數(shù)掃描本身代碼量不大但勝在能暴露模型的兩個隱藏問題RE 閾值定太緊時模型變不可行閾值定太松時 Q 沒有明顯改善這兩個現(xiàn)象直接告訴你模型的可行域有多“窄”數(shù)據(jù)積累得多了判斷資源瓶頸在哪、需求哪里過高比單次最優(yōu)解可靠得多。希望幫到你。本文還有配套的精品資源點擊獲取