劃實(shí)戰(zhàn)：從建模到靈敏度分析，掌握l(shuí)inprog核心應(yīng)用)
1. 項(xiàng)目概述從練習(xí)題到實(shí)戰(zhàn)思維的跨越很多朋友一看到“Matlab 線性規(guī)劃練習(xí)題”可能第一反應(yīng)就是翻開(kāi)教材對(duì)著例題敲幾行代碼得到答案就完事了。我剛開(kāi)始接觸運(yùn)籌學(xué)和Matlab的時(shí)候也是這么干的但后來(lái)在實(shí)際項(xiàng)目中吃了虧才發(fā)現(xiàn)這種“做題”思維和解決真實(shí)問(wèn)題的“建?！彼季S中間隔著一道巨大的鴻溝。這道鴻溝不是靠解出正確答案就能填平的它關(guān)乎你如何定義問(wèn)題、如何處理數(shù)據(jù)、如何解讀結(jié)果以及最關(guān)鍵的一一當(dāng)模型與現(xiàn)實(shí)出現(xiàn)偏差時(shí)你該如何調(diào)整。所以今天我想聊的遠(yuǎn)不止是幾道練習(xí)題的標(biāo)準(zhǔn)解法。我想和你分享的是如何利用Matlab這個(gè)強(qiáng)大的工具把線性規(guī)劃從課本上的數(shù)學(xué)公式變成你手中解決實(shí)際工程、經(jīng)濟(jì)、管理問(wèn)題的“瑞士軍刀”。我們會(huì)從最基礎(chǔ)的模型構(gòu)建講起一步步深入到靈敏度分析、大規(guī)模問(wèn)題處理和那些教科書(shū)里不會(huì)寫(xiě)的“踩坑”經(jīng)驗(yàn)。無(wú)論你是正在備考的學(xué)生還是希望將線性規(guī)劃應(yīng)用于實(shí)際工作的工程師、分析師我相信這些從一線實(shí)踐中總結(jié)出來(lái)的思路和技巧都能讓你少走很多彎路。線性規(guī)劃的核心很簡(jiǎn)單在滿足一組線性等式或不等式約束的條件下最大化或最小化一個(gè)線性目標(biāo)函數(shù)。但它的魅力就在于這個(gè)簡(jiǎn)單的框架可以描述生產(chǎn)計(jì)劃、資源分配、投資組合、運(yùn)輸調(diào)度等成千上萬(wàn)的現(xiàn)實(shí)問(wèn)題。Matlab的linprog函數(shù)就像是一個(gè)封裝好的“求解器黑箱”我們的核心工作其實(shí)是在“黑箱”之外——如何把一團(tuán)亂麻的現(xiàn)實(shí)精準(zhǔn)地翻譯成linprog能聽(tīng)懂的“語(yǔ)言”即標(biāo)準(zhǔn)型并理解它返回的每一個(gè)數(shù)字背后的業(yè)務(wù)意義。2. 線性規(guī)劃模型構(gòu)建與Matlab標(biāo)準(zhǔn)化全解析2.1 理解線性規(guī)劃的標(biāo)準(zhǔn)形式一切求解的起點(diǎn)Matlab的linprog函數(shù)只認(rèn)一種“方言”那就是線性規(guī)劃的標(biāo)準(zhǔn)形式最小值形式。很多初學(xué)者栽的第一個(gè)跟頭就是憑感覺(jué)把問(wèn)題輸進(jìn)去結(jié)果報(bào)錯(cuò)或者得出一個(gè)莫名其妙的答案。我們必須像強(qiáng)迫癥一樣把任何問(wèn)題都先轉(zhuǎn)化為這個(gè)標(biāo)準(zhǔn)形式標(biāo)準(zhǔn)形式最小化 ( f^T x ) 滿足 ( A \cdot x \leq b ) ( Aeq \cdot x beq ) ( lb \leq x \leq ub )其中( x )是我們的決策變量向量。這里有幾個(gè)關(guān)鍵點(diǎn)必須死記硬背目標(biāo)必須是“最小化”如果你的問(wèn)題是最大化利潤(rùn)比如max z 3x1 5x2那么你需要將其轉(zhuǎn)化為min z -3x1 -5x2。linprog最終輸出的目標(biāo)函數(shù)值fval也是這個(gè)最小值所以你需要手動(dòng)對(duì)其取負(fù)號(hào)才能得到原問(wèn)題的最大值。不等式約束默認(rèn)為“≤”矩陣A和向量b代表的是 ( A*x \leq b )。如果你的約束是“≥”比如 ( 2x1 x2 \geq 10 )你需要兩邊乘以 -1轉(zhuǎn)化為 ( -2x1 - x2 \leq -10 )。決策變量默認(rèn)非負(fù)在標(biāo)準(zhǔn)形式中隱含了 ( x \geq 0 ) 的條件這由下界lb來(lái)指定。如果你的變量可以取負(fù)值必須顯式地將對(duì)應(yīng)變量的下界lb設(shè)置為-Inf。注意這是最容易出錯(cuò)的地方之一。我見(jiàn)過(guò)有人把“≥”約束直接寫(xiě)成A*x b傳入linprog結(jié)果求解器把它當(dāng)作“≤”來(lái)處理導(dǎo)致完全錯(cuò)誤的可行域。養(yǎng)成第一步先手工轉(zhuǎn)化為標(biāo)準(zhǔn)型的習(xí)慣能節(jié)省大量調(diào)試時(shí)間。2.2 從文字描述到數(shù)學(xué)模型的實(shí)戰(zhàn)拆解看一個(gè)練習(xí)題“某工廠生產(chǎn)A、B兩種產(chǎn)品生產(chǎn)每件A產(chǎn)品需耗材2kg、工時(shí)1小時(shí)利潤(rùn)3元生產(chǎn)每件B產(chǎn)品需耗材1kg、工時(shí)2小時(shí)利潤(rùn)5元。現(xiàn)有材料100kg工時(shí)80小時(shí)。問(wèn)如何安排生產(chǎn)使利潤(rùn)最大”第一步定義決策變量這看似簡(jiǎn)單卻至關(guān)重要。設(shè) ( x_1 ) 為產(chǎn)品A的產(chǎn)量( x_2 ) 為產(chǎn)品B的產(chǎn)量。變量名最好有明確意義如xA,xB在復(fù)雜問(wèn)題中能極大提升代碼可讀性。第二步構(gòu)建目標(biāo)函數(shù)目標(biāo)是最大化總利潤(rùn)( Max , Z 3x_1 5x_2 )。為了適配linprog我們將其改寫(xiě)為( Min , Z -3x_1 - 5x_2 )。所以Matlab中的目標(biāo)函數(shù)系數(shù)向量f [-3; -5]。第三步提煉約束條件材料約束( 2x_1 x_2 \leq 100 )工時(shí)約束( x_1 2x_2 \leq 80 )非負(fù)約束( x_1 \geq 0, x_2 \geq 0 )這里約束都是“≤”且變量非負(fù)符合標(biāo)準(zhǔn)形式。因此我們可以直接寫(xiě)出不等式約束矩陣A [2, 1; 1, 2]不等式約束右側(cè)向量b [100; 80]下界lb [0; 0](上界ub默認(rèn)為正無(wú)窮無(wú)需指定)第四步識(shí)別等式約束本題沒(méi)有等式約束所以Aeq和beq為空矩陣[]。通過(guò)這個(gè)步驟我們就把一段文字描述嚴(yán)謹(jǐn)?shù)胤g成了linprog(f, A, b, Aeq, beq, lb, ub)函數(shù)所需要的各個(gè)輸入?yún)?shù)。這個(gè)過(guò)程就是數(shù)學(xué)建模的核心。2.3 輸入?yún)?shù)的組織技巧與常見(jiàn)陷阱在Matlab中組織這些矩陣和向量時(shí)有一些細(xì)節(jié)技巧% 推薦的做法清晰、易核對(duì) f [-3; -5]; % 目標(biāo)函數(shù)系數(shù)列向量 A [2, 1; % 第一行材料消耗系數(shù) 1, 2]; % 第二行工時(shí)消耗系數(shù) b [100; 80]; Aeq []; % 無(wú)等式約束顯式賦空 beq []; lb [0; 0]; % 非負(fù)約束 % 調(diào)用linprog求解 [x_opt, fval_opt] linprog(f, A, b, Aeq, beq, lb); % 注意linprog返回的是最小值所以最大利潤(rùn)需要取負(fù) max_profit -fval_opt; fprintf(最優(yōu)生產(chǎn)計(jì)劃生產(chǎn)A產(chǎn)品 %.2f 件生產(chǎn)B產(chǎn)品 %.2f 件。\n, x_opt(1), x_opt(2)); fprintf(最大利潤(rùn)為%.2f 元。\n, max_profit);常見(jiàn)陷阱維度不匹配確保f的長(zhǎng)度決策變量個(gè)數(shù)與A的列數(shù)、Aeq的列數(shù)、lb的長(zhǎng)度完全一致。A的行數(shù)等于不等式約束的個(gè)數(shù)必須與b的長(zhǎng)度一致。這是最常出現(xiàn)的錯(cuò)誤之一。稀疏矩陣問(wèn)題對(duì)于大規(guī)模線性規(guī)劃變量或約束成千上萬(wàn)如果A或Aeq中大部分元素是0使用稀疏矩陣sparse存儲(chǔ)可以極大節(jié)省內(nèi)存和提高求解速度。linprog支持稀疏矩陣輸入。選項(xiàng)配置默認(rèn)的求解器和設(shè)置可能對(duì)某些病態(tài)問(wèn)題不理想。可以通過(guò)optimoptions來(lái)調(diào)整例如顯示迭代過(guò)程、設(shè)置容差或選擇算法。options optimoptions(linprog, Display, iter, Algorithm, dual-simplex); [x_opt, fval_opt] linprog(f, A, b, Aeq, beq, lb, ub, options);‘dual-simplex’對(duì)偶單純形法通常對(duì)重新求解或靈敏度分析更高效。3. 練習(xí)題深度求解與結(jié)果解讀3.1 基礎(chǔ)練習(xí)題求解與可視化分析我們求解上述生產(chǎn)計(jì)劃問(wèn)題。運(yùn)行代碼后得到最優(yōu)解x1 40, x2 20 最優(yōu)目標(biāo)函數(shù)值最小值fval -220 因此最大利潤(rùn)為 220 元。但這只是數(shù)字。為了真正理解我們可以可視化可行域和最優(yōu)解。對(duì)于二維問(wèn)題這是極佳的學(xué)習(xí)手段。% 繪制約束條件 figure; hold on; grid on; % 約束1: 2x1 x2 100 - x2 100 - 2*x1 x1_line linspace(0, 50, 100); x2_line1 100 - 2*x1_line; plot(x1_line, x2_line1, b-, LineWidth, 2); % 約束2: x1 2x2 80 - x2 (80 - x1)/2 x2_line2 (80 - x1_line)/2; plot(x1_line, x2_line2, r-, LineWidth, 2); % 非負(fù)約束即坐標(biāo)軸 xlim([0, 60]); ylim([0, 60]); xlabel(產(chǎn)品A產(chǎn)量 x1); ylabel(產(chǎn)品B產(chǎn)量 x2); % 填充可行域多邊形頂點(diǎn) % 可行域頂點(diǎn)由約束方程兩兩相交及與坐標(biāo)軸相交得到 vertices [0,0; 50,0; 40,20; 0,40]; % 按順序連接構(gòu)成可行域多邊形 fill(vertices(:,1), vertices(:,2), g, FaceAlpha, 0.2); % 繪制目標(biāo)函數(shù)等值線及最優(yōu)解 % 目標(biāo)函數(shù): Z 3x1 5x2 - x2 (Z - 3x1)/5 Z_values [100, 220, 300]; % 繪制幾條等利潤(rùn)線 for Z Z_values x2_obj (Z - 3*x1_line)/5; plot(x1_line, x2_obj, k--); text(x1_line(end), x2_obj(end), sprintf(Z%d, Z)); end % 標(biāo)記最優(yōu)解點(diǎn) plot(x_opt(1), x_opt(2), ro, MarkerSize, 10, MarkerFaceColor, r); text(x_opt(1)2, x_opt(2)2, sprintf(最優(yōu)解 (%.0f,%.0f), x_opt), FontWeight, bold); legend(材料約束, 工時(shí)約束, 可行域, 等利潤(rùn)線, 最優(yōu)解); title(線性規(guī)劃問(wèn)題可行域與最優(yōu)解可視化);通過(guò)這張圖你可以清晰地看到由藍(lán)色和紅色約束線、坐標(biāo)軸圍成的綠色區(qū)域就是可行域所有可能的生產(chǎn)計(jì)劃都落在這個(gè)區(qū)域內(nèi)。黑色虛線是等利潤(rùn)線同一條線上的點(diǎn)利潤(rùn)相同。最優(yōu)解紅點(diǎn)位于可行域的一個(gè)頂點(diǎn)上并且是等利潤(rùn)線在向右上方平移時(shí)利潤(rùn)增加最后離開(kāi)可行域的那個(gè)點(diǎn)。這直觀驗(yàn)證了線性規(guī)劃的一個(gè)核心定理最優(yōu)解如果存在必定在可行域的某個(gè)頂點(diǎn)取得。3.2 結(jié)果的經(jīng)濟(jì)學(xué)與管理學(xué)解讀得到x140, x220和max Z220后工作只完成了一半。更重要的是解讀松弛變量我們可以計(jì)算每個(gè)約束的實(shí)際使用情況。材料2*40 1*20 100kg恰好用完松弛變量為0。工時(shí)1*40 2*20 80小時(shí)也恰好用完松弛變量為0。管理啟示這個(gè)解被稱(chēng)為“緊致的”或“活躍的”因?yàn)樗匈Y源都耗盡了。這意味著材料和工時(shí)都是瓶頸資源任何一項(xiàng)的增加都可能帶來(lái)利潤(rùn)的提升。這自然引出了下一個(gè)關(guān)鍵問(wèn)題如果我能多獲得1kg材料或多1小時(shí)工時(shí)利潤(rùn)能增加多少這就是靈敏度分析影子價(jià)格要回答的問(wèn)題。linprog函數(shù)可以通過(guò)輸出更多參數(shù)來(lái)直接獲得這些信息[x_opt, fval_opt, exitflag, output, lambda] linprog(f, A, b, Aeq, beq, lb); % lambda.ineqlin 給出了不等式約束的影子價(jià)格拉格朗日乘子 % lambda.lower 給出了下界約束的影子價(jià)格 shadow_price_material lambda.ineqlin(1); shadow_price_labor lambda.ineqlin(2); fprintf(材料的影子價(jià)格邊際利潤(rùn)為%.4f 元/千克\n, shadow_price_material); fprintf(工時(shí)的影子價(jià)格邊際利潤(rùn)為%.4f 元/小時(shí)\n, shadow_price_labor);影子價(jià)格是管理層做決策的黃金指標(biāo)。如果市場(chǎng)上材料的采購(gòu)成本低于其影子價(jià)格那么增加材料儲(chǔ)備就是劃算的反之則不應(yīng)增加。3.3 中等難度練習(xí)題混合與比例問(wèn)題來(lái)看一個(gè)更復(fù)雜點(diǎn)的例子“一種飼料由玉米、豆粕、麥麩混合而成。要求每公斤飼料中蛋白質(zhì)至少18%纖維至多8%脂肪至少5%。三種原料的營(yíng)養(yǎng)成分及價(jià)格如下表。求滿足營(yíng)養(yǎng)要求且成本最低的混合方案?！痹系鞍踪|(zhì)(%)纖維(%)脂肪(%)成本(元/kg)玉米8342.0豆粕45614.5麥麩151031.5建模步驟變量設(shè)每公斤飼料中玉米、豆粕、麥麩的含量分別為 ( x_1, x_2, x_3 ) 公斤。注意這里變量代表的是比例或重量總和應(yīng)為1公斤。目標(biāo)最小化成本 ( Min , Z 2.0x_1 4.5x_2 1.5x_3 )。約束營(yíng)養(yǎng)約束注意單位轉(zhuǎn)換例如蛋白質(zhì)含量約束蛋白質(zhì)( 0.08x_1 0.45x_2 0.15x_3 \geq 0.18 ) (至少18%)纖維( 0.03x_1 0.06x_2 0.10x_3 \leq 0.08 ) (至多8%)脂肪( 0.04x_1 0.01x_2 0.03x_3 \geq 0.05 ) (至少5%)比例約束所有原料加起來(lái)正好是一公斤( x_1 x_2 x_3 1 )非負(fù)約束( x_1, x_2, x_3 \geq 0 )轉(zhuǎn)化為Matlab標(biāo)準(zhǔn)型目標(biāo)函數(shù)f [2.0; 4.5; 1.5](已是最小化)不等式約束需要處理“≥”。蛋白質(zhì)和脂肪約束是“≥”需乘以-1。蛋白質(zhì)( -0.08x_1 - 0.45x_2 - 0.15x_3 \leq -0.18 )脂肪( -0.04x_1 - 0.01x_2 - 0.03x_3 \leq -0.05 )纖維( 0.03x_1 0.06x_2 0.10x_3 \leq 0.08 ) 所以A [-0.08, -0.45, -0.15; -0.04, -0.01, -0.03; 0.03, 0.06, 0.10],b [-0.18; -0.05; 0.08]等式約束Aeq [1, 1, 1],beq [1]下界lb [0; 0; 0]求解這個(gè)模型不僅能得到最低成本配方還能分析哪種營(yíng)養(yǎng)要求是“緊約束”以及每種原料的影子價(jià)格從而指導(dǎo)采購(gòu)策略。例如如果豆粕的影子價(jià)格很高說(shuō)明其提供的蛋白質(zhì)非常“珍貴”即使價(jià)格高在配方中也難以被完全替代。4. 高級(jí)應(yīng)用靈敏度分析與參數(shù)規(guī)劃4.1 影子價(jià)格與資源“價(jià)值”評(píng)估靈敏度分析是線性規(guī)劃的靈魂它回答“如果……會(huì)怎樣”的問(wèn)題。linprog返回的lambda結(jié)構(gòu)體就包含了這些信息。lambda.ineqlin: 對(duì)應(yīng)不等式約束 ( A*x \leq b )。其分量的值就是該約束的影子價(jià)格。它表示對(duì)應(yīng)約束的右側(cè)常數(shù) ( b_i ) 每增加一個(gè)單位放松約束最優(yōu)目標(biāo)函數(shù)值這里是成本能改善多少對(duì)于最小化問(wèn)題是減少對(duì)于最大化問(wèn)題是增加。lambda.eqlin: 對(duì)應(yīng)等式約束 ( Aeq*x beq ) 的影子價(jià)格。lambda.lower/upper: 對(duì)應(yīng)變量上下界約束的影子價(jià)格。在前面的飼料例子中求解后查看lambda.ineqlin第一個(gè)值對(duì)應(yīng)蛋白質(zhì)約束可能是正的。這意味著如果蛋白質(zhì)要求從18%略微放寬到17.99%單位飼料成本可以下降多少元。這個(gè)值代表了滿足該營(yíng)養(yǎng)要求的“邊際成本”。如果某個(gè)約束的影子價(jià)格為0通常對(duì)應(yīng)該資源有剩余松弛變量0則說(shuō)明該約束不是當(dāng)前最優(yōu)解的活躍約束稍微調(diào)整它不會(huì)影響最優(yōu)成本。實(shí)操心得影子價(jià)格只在最優(yōu)基不變的一個(gè)小范圍內(nèi)有效這個(gè)范圍就是靈敏度報(bào)告的內(nèi)容。Matlab的linprog沒(méi)有直接提供完整的靈敏度范圍但你可以通過(guò)手動(dòng)微調(diào)b或c的值重新求解來(lái)觀察最優(yōu)解在何時(shí)發(fā)生變化從而估算出這個(gè)有效范圍。對(duì)于關(guān)鍵決策這一步必不可少。4.2 目標(biāo)函數(shù)系數(shù)變化的靈敏度分析除了資源數(shù)量b變化原料價(jià)格目標(biāo)函數(shù)系數(shù)f也可能波動(dòng)。我們需要知道價(jià)格在什么范圍內(nèi)波動(dòng)時(shí)當(dāng)前的最優(yōu)生產(chǎn)方案即哪些產(chǎn)品生產(chǎn)哪些不生產(chǎn)保持不變。對(duì)于簡(jiǎn)單的二維問(wèn)題我們可以通過(guò)圖形觀察等成本線斜率的變化范圍。對(duì)于高維問(wèn)題則需要通過(guò)計(jì)算** Reduced Cost檢驗(yàn)數(shù)** 來(lái)分析。在linprog的求解結(jié)果中非基變量的 Reduced Cost 可以直接反映其經(jīng)濟(jì)意義它表示該變量比如一種當(dāng)前不使用的原料要變得有利可圖其成本需要降低多少對(duì)于最小化問(wèn)題。雖然linprog不直接輸出 Reduced Cost 矩陣但我們可以根據(jù)最優(yōu)解和拉格朗日乘子來(lái)理解對(duì)于一個(gè)取值為0的變量即不在最優(yōu)生產(chǎn)組合中如果其 Reduced Cost 為正最小化問(wèn)題意味著即使該變量的成本再降低一點(diǎn)引入它也不會(huì)讓總成本更低當(dāng)前不用的決策是穩(wěn)健的。4.3 參數(shù)線性規(guī)劃實(shí)戰(zhàn)應(yīng)對(duì)市場(chǎng)波動(dòng)現(xiàn)實(shí)世界中價(jià)格和資源量都不是固定不變的。參數(shù)規(guī)劃研究當(dāng)這些參數(shù)連續(xù)變化時(shí)最優(yōu)解如何變化。我們可以通過(guò)循環(huán)或更高級(jí)的算法來(lái)實(shí)現(xiàn)。例如在飼料問(wèn)題中假設(shè)豆粕價(jià)格 ( c_2 ) 是一個(gè)波動(dòng)參數(shù)。我們想畫(huà)出“最低飼料成本”隨“豆粕價(jià)格”變化的函數(shù)圖像。% 參數(shù)規(guī)劃示例豆粕價(jià)格變化對(duì)總成本的影響 c2_range 3.0:0.1:6.0; % 假設(shè)豆粕價(jià)格在3元到6元之間波動(dòng) optimal_cost zeros(size(c2_range)); optimal_x2 zeros(size(c2_range)); for i 1:length(c2_range) f_param [2.0; c2_range(i); 1.5]; % 更新目標(biāo)函數(shù)系數(shù) [x, fval] linprog(f_param, A, b, Aeq, beq, lb); optimal_cost(i) fval; optimal_x2(i) x(2); % 記錄豆粕的最優(yōu)使用量 end figure; subplot(2,1,1); plot(c2_range, optimal_cost, b-o, LineWidth, 1.5); xlabel(豆粕價(jià)格 (元/kg)); ylabel(最低飼料成本 (元/kg)); title(最低成本隨豆粕價(jià)格變化曲線); grid on; subplot(2,1,2); plot(c2_range, optimal_x2, r-s, LineWidth, 1.5); xlabel(豆粕價(jià)格 (元/kg)); ylabel(最優(yōu)配方中豆粕含量 (kg)); title(豆粕最優(yōu)用量隨其價(jià)格變化曲線); grid on;從圖中你可以清晰地看到成本曲線是一條分段線性凸函數(shù)。在豆粕價(jià)格較低時(shí)成本隨價(jià)格上升緩慢線性增加同時(shí)豆粕用量保持高位。當(dāng)價(jià)格超過(guò)某個(gè)臨界點(diǎn)豆粕用量會(huì)突然下降被其他原料替代成本曲線的斜率也會(huì)發(fā)生變化。這個(gè)臨界點(diǎn)就是最優(yōu)基變更的點(diǎn)。這種分析對(duì)于采購(gòu)策略和風(fēng)險(xiǎn)管理至關(guān)重要。它能告訴你在什么價(jià)格區(qū)間內(nèi)可以維持現(xiàn)有的配方結(jié)構(gòu)當(dāng)價(jià)格突破什么閾值時(shí)就必須考慮調(diào)整配方。5. 大規(guī)模問(wèn)題求解與性能優(yōu)化技巧5.1 稀疏矩陣與模型高效構(gòu)建當(dāng)問(wèn)題規(guī)模變大例如涉及數(shù)百個(gè)倉(cāng)庫(kù)向數(shù)千個(gè)門(mén)店的運(yùn)輸問(wèn)題約束矩陣A會(huì)非常龐大但其中絕大多數(shù)元素是0每個(gè)約束通常只涉及少數(shù)變量。使用全矩陣存儲(chǔ)會(huì)消耗巨大內(nèi)存。Matlab的稀疏矩陣只存儲(chǔ)非零元素的位置和值。構(gòu)建大規(guī)模LP模型時(shí)應(yīng)優(yōu)先采用稀疏格式。% 假設(shè)有一個(gè)1000個(gè)變量500個(gè)約束的問(wèn)題每個(gè)約束平均涉及5個(gè)變量 n_vars 1000; n_ineq 500; % 預(yù)先估算非零元個(gè)數(shù) nnz_estimated n_ineq * 5; % 使用稀疏矩陣構(gòu)造函數(shù)效率遠(yuǎn)高于先創(chuàng)建全矩陣再轉(zhuǎn)換 A_sparse sparse(n_ineq, n_vars); % 先創(chuàng)建全零稀疏矩陣 % 或者更高效地直接使用三元組格式構(gòu)建 rows []; cols []; vals []; % ... 在循環(huán)中填充 rows, cols, vals ... A_sparse sparse(rows, cols, vals, n_ineq, n_vars); f randn(n_vars, 1); % 目標(biāo)函數(shù)系數(shù) b rand(n_ineq, 1) * 100; lb zeros(n_vars, 1); % 調(diào)用linprog傳入稀疏矩陣 options optimoptions(linprog, Display, final, Algorithm, dual-simplex); [x, fval] linprog(f, A_sparse, b, [], [], lb, [], [], options);使用稀疏矩陣后內(nèi)存占用和求解時(shí)間通常會(huì)大幅下降。5.2 算法選擇與求解器選項(xiàng)調(diào)優(yōu)linprog內(nèi)置了多種算法通過(guò)Algorithm選項(xiàng)指定‘dual-simplex’(默認(rèn))對(duì)偶單純形法。對(duì)于重新求解、熱啟動(dòng)、以及約束遠(yuǎn)多于變量的問(wèn)題通常很高效。在進(jìn)行靈敏度分析或參數(shù)規(guī)劃時(shí)如果每次只改變少量數(shù)據(jù)此算法可能更快?！甶nterior-point’內(nèi)點(diǎn)法。對(duì)于大規(guī)模、稠密的問(wèn)題通常比單純形法更快尤其是變量和約束數(shù)量都很大的情況。它總是通過(guò)內(nèi)點(diǎn)路徑逼近最優(yōu)解而單純形法在頂點(diǎn)間移動(dòng)。‘interior-point-legacy’舊版本的內(nèi)點(diǎn)法。選擇建議對(duì)于中小規(guī)模問(wèn)題變量/約束數(shù)10000兩種算法差異不大可以都試試。對(duì)于大規(guī)模問(wèn)題尤其是結(jié)構(gòu)稀疏的優(yōu)先嘗試‘interior-point’。如果你需要頻繁地修改模型參數(shù)如做參數(shù)分析并且每次修改很小使用‘dual-simplex’并配合熱啟動(dòng)可能更優(yōu)。熱啟動(dòng)可以通過(guò)設(shè)置初始點(diǎn)x0來(lái)實(shí)現(xiàn)雖然linprog對(duì)熱啟動(dòng)的支持不如專(zhuān)業(yè)求解器但提供一個(gè)接近最優(yōu)解的初始點(diǎn)有時(shí)能加快收斂。% 使用內(nèi)點(diǎn)法求解大規(guī)模問(wèn)題 options optimoptions(linprog, ... Algorithm, interior-point, ... Display, iter, ... % 顯示迭代信息 OptimalityTolerance, 1e-8, ... % 優(yōu)化容差 ConstraintTolerance, 1e-8); % 約束容差 [x, fval] linprog(f, A, b, Aeq, beq, lb, ub, [], options);5.3 模型調(diào)試與錯(cuò)誤排查實(shí)戰(zhàn)指南即使經(jīng)驗(yàn)豐富建模錯(cuò)誤也時(shí)有發(fā)生。以下是一些快速排查清單Linprog stopped because no point satisfies the constraints.(問(wèn)題不可行)原因約束條件互相矛盾沒(méi)有同時(shí)滿足所有條件的解。排查檢查是否有明顯的矛盾如x 5和x 10。檢查所有“≥”約束是否已正確轉(zhuǎn)換為“≤”。逐步注釋掉部分約束找到導(dǎo)致不可行的那個(gè)約束。使用linprog的‘diagnostics’選項(xiàng)獲取更詳細(xì)的信息。options optimoptions(linprog, Display, iter-detailed);Linprog stopped because the objective function is unbounded below.(問(wèn)題無(wú)界)原因在可行域內(nèi)目標(biāo)函數(shù)值可以無(wú)限減小對(duì)于最小化問(wèn)題。通常是因?yàn)橥浱砑颖匾募s束比如變量沒(méi)有上界且其在目標(biāo)函數(shù)中的系數(shù)為負(fù)。排查檢查是否所有變量都有合理的上下界。檢查目標(biāo)函數(shù)系數(shù)的符號(hào)。在最大化問(wèn)題中忘記對(duì)目標(biāo)函數(shù)取負(fù)號(hào)也可能導(dǎo)致此錯(cuò)誤因?yàn)閘inprog在嘗試最小化一個(gè)正系數(shù)函數(shù)且變量無(wú)上界。得到非整數(shù)解但期望整數(shù)解原因linprog求解的是連續(xù)線性規(guī)劃。如果你需要整數(shù)解如生產(chǎn)數(shù)量不能是小數(shù)那么你面對(duì)的是一個(gè)整數(shù)線性規(guī)劃問(wèn)題需要使用intlinprog函數(shù)。解決將linprog替換為intlinprog并指定哪些變量需要取整。% 假設(shè) x1 和 x2 需要是整數(shù) intcon [1, 2]; % 整數(shù)變量的索引 [x_opt, fval_opt] intlinprog(f, intcon, A, b, Aeq, beq, lb, ub);注意整數(shù)規(guī)劃求解難度遠(yuǎn)大于線性規(guī)劃求解時(shí)間可能很長(zhǎng)甚至對(duì)大規(guī)模問(wèn)題無(wú)法在可接受時(shí)間內(nèi)找到最優(yōu)解。求解速度慢嘗試使用稀疏矩陣。切換算法在‘dual-simplex’和‘interior-point’之間切換。調(diào)整優(yōu)化選項(xiàng)適當(dāng)放寬OptimalityTolerance和ConstraintTolerance如從1e-8調(diào)到1e-6可以加速收斂尤其對(duì)于內(nèi)點(diǎn)法。檢查模型規(guī)模如果可能?chē)L試從業(yè)務(wù)邏輯上簡(jiǎn)化模型減少不必要的變量或約束。6. 從練習(xí)題到實(shí)際項(xiàng)目思維轉(zhuǎn)變與經(jīng)驗(yàn)分享做了這么多練習(xí)題最終目的是為了應(yīng)用于實(shí)際。在實(shí)際項(xiàng)目中線性規(guī)劃的應(yīng)用遠(yuǎn)比課本習(xí)題復(fù)雜。首先數(shù)據(jù)準(zhǔn)備就是一大挑戰(zhàn)。練習(xí)題里的數(shù)據(jù)都是干凈的、給定的?，F(xiàn)實(shí)中原料價(jià)格、資源消耗系數(shù)、市場(chǎng)需求等數(shù)據(jù)都來(lái)自數(shù)據(jù)庫(kù)、報(bào)表或預(yù)測(cè)模型可能存在缺失、異常或不確定性。在構(gòu)建模型前必須進(jìn)行嚴(yán)格的數(shù)據(jù)清洗和驗(yàn)證。例如單位是否統(tǒng)一歷史消耗系數(shù)是否適用于未來(lái)我曾在做一個(gè)生產(chǎn)排程項(xiàng)目時(shí)因?yàn)檎`將“噸”和“公斤”的數(shù)據(jù)混用導(dǎo)致整個(gè)模型失效教訓(xùn)深刻。其次模型的維護(hù)和更新至關(guān)重要。一個(gè)線性規(guī)劃模型不是一勞永逸的。市場(chǎng)條件、生產(chǎn)工藝、公司政策都在變模型必須能方便地更新參數(shù)甚至調(diào)整結(jié)構(gòu)。這就要求我們?cè)诰幊虝r(shí)要有良好的習(xí)慣將數(shù)據(jù)f,A,b等與模型邏輯分離最好從外部文件如Excel、CSV讀取數(shù)據(jù)。這樣業(yè)務(wù)人員只需更新數(shù)據(jù)文件而無(wú)需觸碰代碼。% 良好的實(shí)踐從外部文件讀取數(shù)據(jù) data readmatrix(production_data.csv); f data(1, :); % 第一行是目標(biāo)系數(shù) A data(2:end-1, :); % 中間行是不等式約束矩陣 b data(end, :); % 最后一行是約束右側(cè)值 % ... 然后調(diào)用linprog最后結(jié)果的呈現(xiàn)和解釋需要技巧。你不能只給老板看一堆x1123.456這樣的數(shù)字。你需要將其轉(zhuǎn)化為業(yè)務(wù)語(yǔ)言“下個(gè)月建議將A生產(chǎn)線產(chǎn)能的60%用于生產(chǎn)產(chǎn)品X這樣可以預(yù)期提升毛利15%但需要提前儲(chǔ)備B原料20噸。” 同時(shí)一定要附上關(guān)鍵的前提假設(shè)和靈敏度分析結(jié)論“這個(gè)建議是基于產(chǎn)品Y價(jià)格穩(wěn)定在每單位100元的預(yù)測(cè)。如果其價(jià)格下跌超過(guò)5%方案需要重新評(píng)估?！本€性規(guī)劃是一個(gè)強(qiáng)大的工具但它只是一個(gè)工具。它的輸出質(zhì)量完全取決于輸入模型的質(zhì)量和對(duì)現(xiàn)實(shí)問(wèn)題的抽象能力。練習(xí)題教會(huì)我們?nèi)绾问褂眠@個(gè)工具而真正的藝術(shù)在于如何用它來(lái)刻畫(huà)和解決那些模糊、復(fù)雜、動(dòng)態(tài)的真實(shí)世界問(wèn)題。這需要不斷的練習(xí)、反思以及最重要的——與實(shí)際問(wèn)題深度碰撞。 Matlab線性規(guī)劃實(shí)戰(zhàn)：從建模到靈敏度分析，掌握l(shuí)inprog核心應(yīng)用
1. 項(xiàng)目概述從練習(xí)題到實(shí)戰(zhàn)思維的跨越很多朋友一看到“Matlab 線性規(guī)劃練習(xí)題”可能第一反應(yīng)就是翻開(kāi)教材對(duì)著例題敲幾行代碼得到答案就完事了。我剛開(kāi)始接觸運(yùn)籌學(xué)和Matlab的時(shí)候也是這么干的但后來(lái)在實(shí)際項(xiàng)目中吃了虧才發(fā)現(xiàn)這種“做題”思維和解決真實(shí)問(wèn)題的“建?！彼季S中間隔著一道巨大的鴻溝。這道鴻溝不是靠解出正確答案就能填平的它關(guān)乎你如何定義問(wèn)題、如何處理數(shù)據(jù)、如何解讀結(jié)果以及最關(guān)鍵的一一當(dāng)模型與現(xiàn)實(shí)出現(xiàn)偏差時(shí)你該如何調(diào)整。所以今天我想聊的遠(yuǎn)不止是幾道練習(xí)題的標(biāo)準(zhǔn)解法。我想和你分享的是如何利用Matlab這個(gè)強(qiáng)大的工具把線性規(guī)劃從課本上的數(shù)學(xué)公式變成你手中解決實(shí)際工程、經(jīng)濟(jì)、管理問(wèn)題的“瑞士軍刀”。我們會(huì)從最基礎(chǔ)的模型構(gòu)建講起一步步深入到靈敏度分析、大規(guī)模問(wèn)題處理和那些教科書(shū)里不會(huì)寫(xiě)的“踩坑”經(jīng)驗(yàn)。無(wú)論你是正在備考的學(xué)生還是希望將線性規(guī)劃應(yīng)用于實(shí)際工作的工程師、分析師我相信這些從一線實(shí)踐中總結(jié)出來(lái)的思路和技巧都能讓你少走很多彎路。線性規(guī)劃的核心很簡(jiǎn)單在滿足一組線性等式或不等式約束的條件下最大化或最小化一個(gè)線性目標(biāo)函數(shù)。但它的魅力就在于這個(gè)簡(jiǎn)單的框架可以描述生產(chǎn)計(jì)劃、資源分配、投資組合、運(yùn)輸調(diào)度等成千上萬(wàn)的現(xiàn)實(shí)問(wèn)題。Matlab的linprog函數(shù)就像是一個(gè)封裝好的“求解器黑箱”我們的核心工作其實(shí)是在“黑箱”之外——如何把一團(tuán)亂麻的現(xiàn)實(shí)精準(zhǔn)地翻譯成linprog能聽(tīng)懂的“語(yǔ)言”即標(biāo)準(zhǔn)型并理解它返回的每一個(gè)數(shù)字背后的業(yè)務(wù)意義。2. 線性規(guī)劃模型構(gòu)建與Matlab標(biāo)準(zhǔn)化全解析2.1 理解線性規(guī)劃的標(biāo)準(zhǔn)形式一切求解的起點(diǎn)Matlab的linprog函數(shù)只認(rèn)一種“方言”那就是線性規(guī)劃的標(biāo)準(zhǔn)形式最小值形式。很多初學(xué)者栽的第一個(gè)跟頭就是憑感覺(jué)把問(wèn)題輸進(jìn)去結(jié)果報(bào)錯(cuò)或者得出一個(gè)莫名其妙的答案。我們必須像強(qiáng)迫癥一樣把任何問(wèn)題都先轉(zhuǎn)化為這個(gè)標(biāo)準(zhǔn)形式標(biāo)準(zhǔn)形式最小化 ( f^T x ) 滿足 ( A \cdot x \leq b ) ( Aeq \cdot x beq ) ( lb \leq x \leq ub )其中( x )是我們的決策變量向量。這里有幾個(gè)關(guān)鍵點(diǎn)必須死記硬背目標(biāo)必須是“最小化”如果你的問(wèn)題是最大化利潤(rùn)比如max z 3x1 5x2那么你需要將其轉(zhuǎn)化為min z -3x1 -5x2。linprog最終輸出的目標(biāo)函數(shù)值fval也是這個(gè)最小值所以你需要手動(dòng)對(duì)其取負(fù)號(hào)才能得到原問(wèn)題的最大值。不等式約束默認(rèn)為“≤”矩陣A和向量b代表的是 ( A*x \leq b )。如果你的約束是“≥”比如 ( 2x1 x2 \geq 10 )你需要兩邊乘以 -1轉(zhuǎn)化為 ( -2x1 - x2 \leq -10 )。決策變量默認(rèn)非負(fù)在標(biāo)準(zhǔn)形式中隱含了 ( x \geq 0 ) 的條件這由下界lb來(lái)指定。如果你的變量可以取負(fù)值必須顯式地將對(duì)應(yīng)變量的下界lb設(shè)置為-Inf。注意這是最容易出錯(cuò)的地方之一。我見(jiàn)過(guò)有人把“≥”約束直接寫(xiě)成A*x b傳入linprog結(jié)果求解器把它當(dāng)作“≤”來(lái)處理導(dǎo)致完全錯(cuò)誤的可行域。養(yǎng)成第一步先手工轉(zhuǎn)化為標(biāo)準(zhǔn)型的習(xí)慣能節(jié)省大量調(diào)試時(shí)間。2.2 從文字描述到數(shù)學(xué)模型的實(shí)戰(zhàn)拆解看一個(gè)練習(xí)題“某工廠生產(chǎn)A、B兩種產(chǎn)品生產(chǎn)每件A產(chǎn)品需耗材2kg、工時(shí)1小時(shí)利潤(rùn)3元生產(chǎn)每件B產(chǎn)品需耗材1kg、工時(shí)2小時(shí)利潤(rùn)5元。現(xiàn)有材料100kg工時(shí)80小時(shí)。問(wèn)如何安排生產(chǎn)使利潤(rùn)最大”第一步定義決策變量這看似簡(jiǎn)單卻至關(guān)重要。設(shè) ( x_1 ) 為產(chǎn)品A的產(chǎn)量( x_2 ) 為產(chǎn)品B的產(chǎn)量。變量名最好有明確意義如xA,xB在復(fù)雜問(wèn)題中能極大提升代碼可讀性。第二步構(gòu)建目標(biāo)函數(shù)目標(biāo)是最大化總利潤(rùn)( Max , Z 3x_1 5x_2 )。為了適配linprog我們將其改寫(xiě)為( Min , Z -3x_1 - 5x_2 )。所以Matlab中的目標(biāo)函數(shù)系數(shù)向量f [-3; -5]。第三步提煉約束條件材料約束( 2x_1 x_2 \leq 100 )工時(shí)約束( x_1 2x_2 \leq 80 )非負(fù)約束( x_1 \geq 0, x_2 \geq 0 )這里約束都是“≤”且變量非負(fù)符合標(biāo)準(zhǔn)形式。因此我們可以直接寫(xiě)出不等式約束矩陣A [2, 1; 1, 2]不等式約束右側(cè)向量b [100; 80]下界lb [0; 0](上界ub默認(rèn)為正無(wú)窮無(wú)需指定)第四步識(shí)別等式約束本題沒(méi)有等式約束所以Aeq和beq為空矩陣[]。通過(guò)這個(gè)步驟我們就把一段文字描述嚴(yán)謹(jǐn)?shù)胤g成了linprog(f, A, b, Aeq, beq, lb, ub)函數(shù)所需要的各個(gè)輸入?yún)?shù)。這個(gè)過(guò)程就是數(shù)學(xué)建模的核心。2.3 輸入?yún)?shù)的組織技巧與常見(jiàn)陷阱在Matlab中組織這些矩陣和向量時(shí)有一些細(xì)節(jié)技巧% 推薦的做法清晰、易核對(duì) f [-3; -5]; % 目標(biāo)函數(shù)系數(shù)列向量 A [2, 1; % 第一行材料消耗系數(shù) 1, 2]; % 第二行工時(shí)消耗系數(shù) b [100; 80]; Aeq []; % 無(wú)等式約束顯式賦空 beq []; lb [0; 0]; % 非負(fù)約束 % 調(diào)用linprog求解 [x_opt, fval_opt] linprog(f, A, b, Aeq, beq, lb); % 注意linprog返回的是最小值所以最大利潤(rùn)需要取負(fù) max_profit -fval_opt; fprintf(最優(yōu)生產(chǎn)計(jì)劃生產(chǎn)A產(chǎn)品 %.2f 件生產(chǎn)B產(chǎn)品 %.2f 件。\n, x_opt(1), x_opt(2)); fprintf(最大利潤(rùn)為%.2f 元。\n, max_profit);常見(jiàn)陷阱維度不匹配確保f的長(zhǎng)度決策變量個(gè)數(shù)與A的列數(shù)、Aeq的列數(shù)、lb的長(zhǎng)度完全一致。A的行數(shù)等于不等式約束的個(gè)數(shù)必須與b的長(zhǎng)度一致。這是最常出現(xiàn)的錯(cuò)誤之一。稀疏矩陣問(wèn)題對(duì)于大規(guī)模線性規(guī)劃變量或約束成千上萬(wàn)如果A或Aeq中大部分元素是0使用稀疏矩陣sparse存儲(chǔ)可以極大節(jié)省內(nèi)存和提高求解速度。linprog支持稀疏矩陣輸入。選項(xiàng)配置默認(rèn)的求解器和設(shè)置可能對(duì)某些病態(tài)問(wèn)題不理想。可以通過(guò)optimoptions來(lái)調(diào)整例如顯示迭代過(guò)程、設(shè)置容差或選擇算法。options optimoptions(linprog, Display, iter, Algorithm, dual-simplex); [x_opt, fval_opt] linprog(f, A, b, Aeq, beq, lb, ub, options);‘dual-simplex’對(duì)偶單純形法通常對(duì)重新求解或靈敏度分析更高效。3. 練習(xí)題深度求解與結(jié)果解讀3.1 基礎(chǔ)練習(xí)題求解與可視化分析我們求解上述生產(chǎn)計(jì)劃問(wèn)題。運(yùn)行代碼后得到最優(yōu)解x1 40, x2 20 最優(yōu)目標(biāo)函數(shù)值最小值fval -220 因此最大利潤(rùn)為 220 元。但這只是數(shù)字。為了真正理解我們可以可視化可行域和最優(yōu)解。對(duì)于二維問(wèn)題這是極佳的學(xué)習(xí)手段。% 繪制約束條件 figure; hold on; grid on; % 約束1: 2x1 x2 100 - x2 100 - 2*x1 x1_line linspace(0, 50, 100); x2_line1 100 - 2*x1_line; plot(x1_line, x2_line1, b-, LineWidth, 2); % 約束2: x1 2x2 80 - x2 (80 - x1)/2 x2_line2 (80 - x1_line)/2; plot(x1_line, x2_line2, r-, LineWidth, 2); % 非負(fù)約束即坐標(biāo)軸 xlim([0, 60]); ylim([0, 60]); xlabel(產(chǎn)品A產(chǎn)量 x1); ylabel(產(chǎn)品B產(chǎn)量 x2); % 填充可行域多邊形頂點(diǎn) % 可行域頂點(diǎn)由約束方程兩兩相交及與坐標(biāo)軸相交得到 vertices [0,0; 50,0; 40,20; 0,40]; % 按順序連接構(gòu)成可行域多邊形 fill(vertices(:,1), vertices(:,2), g, FaceAlpha, 0.2); % 繪制目標(biāo)函數(shù)等值線及最優(yōu)解 % 目標(biāo)函數(shù): Z 3x1 5x2 - x2 (Z - 3x1)/5 Z_values [100, 220, 300]; % 繪制幾條等利潤(rùn)線 for Z Z_values x2_obj (Z - 3*x1_line)/5; plot(x1_line, x2_obj, k--); text(x1_line(end), x2_obj(end), sprintf(Z%d, Z)); end % 標(biāo)記最優(yōu)解點(diǎn) plot(x_opt(1), x_opt(2), ro, MarkerSize, 10, MarkerFaceColor, r); text(x_opt(1)2, x_opt(2)2, sprintf(最優(yōu)解 (%.0f,%.0f), x_opt), FontWeight, bold); legend(材料約束, 工時(shí)約束, 可行域, 等利潤(rùn)線, 最優(yōu)解); title(線性規(guī)劃問(wèn)題可行域與最優(yōu)解可視化);通過(guò)這張圖你可以清晰地看到由藍(lán)色和紅色約束線、坐標(biāo)軸圍成的綠色區(qū)域就是可行域所有可能的生產(chǎn)計(jì)劃都落在這個(gè)區(qū)域內(nèi)。黑色虛線是等利潤(rùn)線同一條線上的點(diǎn)利潤(rùn)相同。最優(yōu)解紅點(diǎn)位于可行域的一個(gè)頂點(diǎn)上并且是等利潤(rùn)線在向右上方平移時(shí)利潤(rùn)增加最后離開(kāi)可行域的那個(gè)點(diǎn)。這直觀驗(yàn)證了線性規(guī)劃的一個(gè)核心定理最優(yōu)解如果存在必定在可行域的某個(gè)頂點(diǎn)取得。3.2 結(jié)果的經(jīng)濟(jì)學(xué)與管理學(xué)解讀得到x140, x220和max Z220后工作只完成了一半。更重要的是解讀松弛變量我們可以計(jì)算每個(gè)約束的實(shí)際使用情況。材料2*40 1*20 100kg恰好用完松弛變量為0。工時(shí)1*40 2*20 80小時(shí)也恰好用完松弛變量為0。管理啟示這個(gè)解被稱(chēng)為“緊致的”或“活躍的”因?yàn)樗匈Y源都耗盡了。這意味著材料和工時(shí)都是瓶頸資源任何一項(xiàng)的增加都可能帶來(lái)利潤(rùn)的提升。這自然引出了下一個(gè)關(guān)鍵問(wèn)題如果我能多獲得1kg材料或多1小時(shí)工時(shí)利潤(rùn)能增加多少這就是靈敏度分析影子價(jià)格要回答的問(wèn)題。linprog函數(shù)可以通過(guò)輸出更多參數(shù)來(lái)直接獲得這些信息[x_opt, fval_opt, exitflag, output, lambda] linprog(f, A, b, Aeq, beq, lb); % lambda.ineqlin 給出了不等式約束的影子價(jià)格拉格朗日乘子 % lambda.lower 給出了下界約束的影子價(jià)格 shadow_price_material lambda.ineqlin(1); shadow_price_labor lambda.ineqlin(2); fprintf(材料的影子價(jià)格邊際利潤(rùn)為%.4f 元/千克\n, shadow_price_material); fprintf(工時(shí)的影子價(jià)格邊際利潤(rùn)為%.4f 元/小時(shí)\n, shadow_price_labor);影子價(jià)格是管理層做決策的黃金指標(biāo)。如果市場(chǎng)上材料的采購(gòu)成本低于其影子價(jià)格那么增加材料儲(chǔ)備就是劃算的反之則不應(yīng)增加。3.3 中等難度練習(xí)題混合與比例問(wèn)題來(lái)看一個(gè)更復(fù)雜點(diǎn)的例子“一種飼料由玉米、豆粕、麥麩混合而成。要求每公斤飼料中蛋白質(zhì)至少18%纖維至多8%脂肪至少5%。三種原料的營(yíng)養(yǎng)成分及價(jià)格如下表。求滿足營(yíng)養(yǎng)要求且成本最低的混合方案?！痹系鞍踪|(zhì)(%)纖維(%)脂肪(%)成本(元/kg)玉米8342.0豆粕45614.5麥麩151031.5建模步驟變量設(shè)每公斤飼料中玉米、豆粕、麥麩的含量分別為 ( x_1, x_2, x_3 ) 公斤。注意這里變量代表的是比例或重量總和應(yīng)為1公斤。目標(biāo)最小化成本 ( Min , Z 2.0x_1 4.5x_2 1.5x_3 )。約束營(yíng)養(yǎng)約束注意單位轉(zhuǎn)換例如蛋白質(zhì)含量約束蛋白質(zhì)( 0.08x_1 0.45x_2 0.15x_3 \geq 0.18 ) (至少18%)纖維( 0.03x_1 0.06x_2 0.10x_3 \leq 0.08 ) (至多8%)脂肪( 0.04x_1 0.01x_2 0.03x_3 \geq 0.05 ) (至少5%)比例約束所有原料加起來(lái)正好是一公斤( x_1 x_2 x_3 1 )非負(fù)約束( x_1, x_2, x_3 \geq 0 )轉(zhuǎn)化為Matlab標(biāo)準(zhǔn)型目標(biāo)函數(shù)f [2.0; 4.5; 1.5](已是最小化)不等式約束需要處理“≥”。蛋白質(zhì)和脂肪約束是“≥”需乘以-1。蛋白質(zhì)( -0.08x_1 - 0.45x_2 - 0.15x_3 \leq -0.18 )脂肪( -0.04x_1 - 0.01x_2 - 0.03x_3 \leq -0.05 )纖維( 0.03x_1 0.06x_2 0.10x_3 \leq 0.08 ) 所以A [-0.08, -0.45, -0.15; -0.04, -0.01, -0.03; 0.03, 0.06, 0.10],b [-0.18; -0.05; 0.08]等式約束Aeq [1, 1, 1],beq [1]下界lb [0; 0; 0]求解這個(gè)模型不僅能得到最低成本配方還能分析哪種營(yíng)養(yǎng)要求是“緊約束”以及每種原料的影子價(jià)格從而指導(dǎo)采購(gòu)策略。例如如果豆粕的影子價(jià)格很高說(shuō)明其提供的蛋白質(zhì)非常“珍貴”即使價(jià)格高在配方中也難以被完全替代。4. 高級(jí)應(yīng)用靈敏度分析與參數(shù)規(guī)劃4.1 影子價(jià)格與資源“價(jià)值”評(píng)估靈敏度分析是線性規(guī)劃的靈魂它回答“如果……會(huì)怎樣”的問(wèn)題。linprog返回的lambda結(jié)構(gòu)體就包含了這些信息。lambda.ineqlin: 對(duì)應(yīng)不等式約束 ( A*x \leq b )。其分量的值就是該約束的影子價(jià)格。它表示對(duì)應(yīng)約束的右側(cè)常數(shù) ( b_i ) 每增加一個(gè)單位放松約束最優(yōu)目標(biāo)函數(shù)值這里是成本能改善多少對(duì)于最小化問(wèn)題是減少對(duì)于最大化問(wèn)題是增加。lambda.eqlin: 對(duì)應(yīng)等式約束 ( Aeq*x beq ) 的影子價(jià)格。lambda.lower/upper: 對(duì)應(yīng)變量上下界約束的影子價(jià)格。在前面的飼料例子中求解后查看lambda.ineqlin第一個(gè)值對(duì)應(yīng)蛋白質(zhì)約束可能是正的。這意味著如果蛋白質(zhì)要求從18%略微放寬到17.99%單位飼料成本可以下降多少元。這個(gè)值代表了滿足該營(yíng)養(yǎng)要求的“邊際成本”。如果某個(gè)約束的影子價(jià)格為0通常對(duì)應(yīng)該資源有剩余松弛變量0則說(shuō)明該約束不是當(dāng)前最優(yōu)解的活躍約束稍微調(diào)整它不會(huì)影響最優(yōu)成本。實(shí)操心得影子價(jià)格只在最優(yōu)基不變的一個(gè)小范圍內(nèi)有效這個(gè)范圍就是靈敏度報(bào)告的內(nèi)容。Matlab的linprog沒(méi)有直接提供完整的靈敏度范圍但你可以通過(guò)手動(dòng)微調(diào)b或c的值重新求解來(lái)觀察最優(yōu)解在何時(shí)發(fā)生變化從而估算出這個(gè)有效范圍。對(duì)于關(guān)鍵決策這一步必不可少。4.2 目標(biāo)函數(shù)系數(shù)變化的靈敏度分析除了資源數(shù)量b變化原料價(jià)格目標(biāo)函數(shù)系數(shù)f也可能波動(dòng)。我們需要知道價(jià)格在什么范圍內(nèi)波動(dòng)時(shí)當(dāng)前的最優(yōu)生產(chǎn)方案即哪些產(chǎn)品生產(chǎn)哪些不生產(chǎn)保持不變。對(duì)于簡(jiǎn)單的二維問(wèn)題我們可以通過(guò)圖形觀察等成本線斜率的變化范圍。對(duì)于高維問(wèn)題則需要通過(guò)計(jì)算** Reduced Cost檢驗(yàn)數(shù)** 來(lái)分析。在linprog的求解結(jié)果中非基變量的 Reduced Cost 可以直接反映其經(jīng)濟(jì)意義它表示該變量比如一種當(dāng)前不使用的原料要變得有利可圖其成本需要降低多少對(duì)于最小化問(wèn)題。雖然linprog不直接輸出 Reduced Cost 矩陣但我們可以根據(jù)最優(yōu)解和拉格朗日乘子來(lái)理解對(duì)于一個(gè)取值為0的變量即不在最優(yōu)生產(chǎn)組合中如果其 Reduced Cost 為正最小化問(wèn)題意味著即使該變量的成本再降低一點(diǎn)引入它也不會(huì)讓總成本更低當(dāng)前不用的決策是穩(wěn)健的。4.3 參數(shù)線性規(guī)劃實(shí)戰(zhàn)應(yīng)對(duì)市場(chǎng)波動(dòng)現(xiàn)實(shí)世界中價(jià)格和資源量都不是固定不變的。參數(shù)規(guī)劃研究當(dāng)這些參數(shù)連續(xù)變化時(shí)最優(yōu)解如何變化。我們可以通過(guò)循環(huán)或更高級(jí)的算法來(lái)實(shí)現(xiàn)。例如在飼料問(wèn)題中假設(shè)豆粕價(jià)格 ( c_2 ) 是一個(gè)波動(dòng)參數(shù)。我們想畫(huà)出“最低飼料成本”隨“豆粕價(jià)格”變化的函數(shù)圖像。% 參數(shù)規(guī)劃示例豆粕價(jià)格變化對(duì)總成本的影響 c2_range 3.0:0.1:6.0; % 假設(shè)豆粕價(jià)格在3元到6元之間波動(dòng) optimal_cost zeros(size(c2_range)); optimal_x2 zeros(size(c2_range)); for i 1:length(c2_range) f_param [2.0; c2_range(i); 1.5]; % 更新目標(biāo)函數(shù)系數(shù) [x, fval] linprog(f_param, A, b, Aeq, beq, lb); optimal_cost(i) fval; optimal_x2(i) x(2); % 記錄豆粕的最優(yōu)使用量 end figure; subplot(2,1,1); plot(c2_range, optimal_cost, b-o, LineWidth, 1.5); xlabel(豆粕價(jià)格 (元/kg)); ylabel(最低飼料成本 (元/kg)); title(最低成本隨豆粕價(jià)格變化曲線); grid on; subplot(2,1,2); plot(c2_range, optimal_x2, r-s, LineWidth, 1.5); xlabel(豆粕價(jià)格 (元/kg)); ylabel(最優(yōu)配方中豆粕含量 (kg)); title(豆粕最優(yōu)用量隨其價(jià)格變化曲線); grid on;從圖中你可以清晰地看到成本曲線是一條分段線性凸函數(shù)。在豆粕價(jià)格較低時(shí)成本隨價(jià)格上升緩慢線性增加同時(shí)豆粕用量保持高位。當(dāng)價(jià)格超過(guò)某個(gè)臨界點(diǎn)豆粕用量會(huì)突然下降被其他原料替代成本曲線的斜率也會(huì)發(fā)生變化。這個(gè)臨界點(diǎn)就是最優(yōu)基變更的點(diǎn)。這種分析對(duì)于采購(gòu)策略和風(fēng)險(xiǎn)管理至關(guān)重要。它能告訴你在什么價(jià)格區(qū)間內(nèi)可以維持現(xiàn)有的配方結(jié)構(gòu)當(dāng)價(jià)格突破什么閾值時(shí)就必須考慮調(diào)整配方。5. 大規(guī)模問(wèn)題求解與性能優(yōu)化技巧5.1 稀疏矩陣與模型高效構(gòu)建當(dāng)問(wèn)題規(guī)模變大例如涉及數(shù)百個(gè)倉(cāng)庫(kù)向數(shù)千個(gè)門(mén)店的運(yùn)輸問(wèn)題約束矩陣A會(huì)非常龐大但其中絕大多數(shù)元素是0每個(gè)約束通常只涉及少數(shù)變量。使用全矩陣存儲(chǔ)會(huì)消耗巨大內(nèi)存。Matlab的稀疏矩陣只存儲(chǔ)非零元素的位置和值。構(gòu)建大規(guī)模LP模型時(shí)應(yīng)優(yōu)先采用稀疏格式。% 假設(shè)有一個(gè)1000個(gè)變量500個(gè)約束的問(wèn)題每個(gè)約束平均涉及5個(gè)變量 n_vars 1000; n_ineq 500; % 預(yù)先估算非零元個(gè)數(shù) nnz_estimated n_ineq * 5; % 使用稀疏矩陣構(gòu)造函數(shù)效率遠(yuǎn)高于先創(chuàng)建全矩陣再轉(zhuǎn)換 A_sparse sparse(n_ineq, n_vars); % 先創(chuàng)建全零稀疏矩陣 % 或者更高效地直接使用三元組格式構(gòu)建 rows []; cols []; vals []; % ... 在循環(huán)中填充 rows, cols, vals ... A_sparse sparse(rows, cols, vals, n_ineq, n_vars); f randn(n_vars, 1); % 目標(biāo)函數(shù)系數(shù) b rand(n_ineq, 1) * 100; lb zeros(n_vars, 1); % 調(diào)用linprog傳入稀疏矩陣 options optimoptions(linprog, Display, final, Algorithm, dual-simplex); [x, fval] linprog(f, A_sparse, b, [], [], lb, [], [], options);使用稀疏矩陣后內(nèi)存占用和求解時(shí)間通常會(huì)大幅下降。5.2 算法選擇與求解器選項(xiàng)調(diào)優(yōu)linprog內(nèi)置了多種算法通過(guò)Algorithm選項(xiàng)指定‘dual-simplex’(默認(rèn))對(duì)偶單純形法。對(duì)于重新求解、熱啟動(dòng)、以及約束遠(yuǎn)多于變量的問(wèn)題通常很高效。在進(jìn)行靈敏度分析或參數(shù)規(guī)劃時(shí)如果每次只改變少量數(shù)據(jù)此算法可能更快?！甶nterior-point’內(nèi)點(diǎn)法。對(duì)于大規(guī)模、稠密的問(wèn)題通常比單純形法更快尤其是變量和約束數(shù)量都很大的情況。它總是通過(guò)內(nèi)點(diǎn)路徑逼近最優(yōu)解而單純形法在頂點(diǎn)間移動(dòng)。‘interior-point-legacy’舊版本的內(nèi)點(diǎn)法。選擇建議對(duì)于中小規(guī)模問(wèn)題變量/約束數(shù)10000兩種算法差異不大可以都試試。對(duì)于大規(guī)模問(wèn)題尤其是結(jié)構(gòu)稀疏的優(yōu)先嘗試‘interior-point’。如果你需要頻繁地修改模型參數(shù)如做參數(shù)分析并且每次修改很小使用‘dual-simplex’并配合熱啟動(dòng)可能更優(yōu)。熱啟動(dòng)可以通過(guò)設(shè)置初始點(diǎn)x0來(lái)實(shí)現(xiàn)雖然linprog對(duì)熱啟動(dòng)的支持不如專(zhuān)業(yè)求解器但提供一個(gè)接近最優(yōu)解的初始點(diǎn)有時(shí)能加快收斂。% 使用內(nèi)點(diǎn)法求解大規(guī)模問(wèn)題 options optimoptions(linprog, ... Algorithm, interior-point, ... Display, iter, ... % 顯示迭代信息 OptimalityTolerance, 1e-8, ... % 優(yōu)化容差 ConstraintTolerance, 1e-8); % 約束容差 [x, fval] linprog(f, A, b, Aeq, beq, lb, ub, [], options);5.3 模型調(diào)試與錯(cuò)誤排查實(shí)戰(zhàn)指南即使經(jīng)驗(yàn)豐富建模錯(cuò)誤也時(shí)有發(fā)生。以下是一些快速排查清單Linprog stopped because no point satisfies the constraints.(問(wèn)題不可行)原因約束條件互相矛盾沒(méi)有同時(shí)滿足所有條件的解。排查檢查是否有明顯的矛盾如x 5和x 10。檢查所有“≥”約束是否已正確轉(zhuǎn)換為“≤”。逐步注釋掉部分約束找到導(dǎo)致不可行的那個(gè)約束。使用linprog的‘diagnostics’選項(xiàng)獲取更詳細(xì)的信息。options optimoptions(linprog, Display, iter-detailed);Linprog stopped because the objective function is unbounded below.(問(wèn)題無(wú)界)原因在可行域內(nèi)目標(biāo)函數(shù)值可以無(wú)限減小對(duì)于最小化問(wèn)題。通常是因?yàn)橥浱砑颖匾募s束比如變量沒(méi)有上界且其在目標(biāo)函數(shù)中的系數(shù)為負(fù)。排查檢查是否所有變量都有合理的上下界。檢查目標(biāo)函數(shù)系數(shù)的符號(hào)。在最大化問(wèn)題中忘記對(duì)目標(biāo)函數(shù)取負(fù)號(hào)也可能導(dǎo)致此錯(cuò)誤因?yàn)閘inprog在嘗試最小化一個(gè)正系數(shù)函數(shù)且變量無(wú)上界。得到非整數(shù)解但期望整數(shù)解原因linprog求解的是連續(xù)線性規(guī)劃。如果你需要整數(shù)解如生產(chǎn)數(shù)量不能是小數(shù)那么你面對(duì)的是一個(gè)整數(shù)線性規(guī)劃問(wèn)題需要使用intlinprog函數(shù)。解決將linprog替換為intlinprog并指定哪些變量需要取整。% 假設(shè) x1 和 x2 需要是整數(shù) intcon [1, 2]; % 整數(shù)變量的索引 [x_opt, fval_opt] intlinprog(f, intcon, A, b, Aeq, beq, lb, ub);注意整數(shù)規(guī)劃求解難度遠(yuǎn)大于線性規(guī)劃求解時(shí)間可能很長(zhǎng)甚至對(duì)大規(guī)模問(wèn)題無(wú)法在可接受時(shí)間內(nèi)找到最優(yōu)解。求解速度慢嘗試使用稀疏矩陣。切換算法在‘dual-simplex’和‘interior-point’之間切換。調(diào)整優(yōu)化選項(xiàng)適當(dāng)放寬OptimalityTolerance和ConstraintTolerance如從1e-8調(diào)到1e-6可以加速收斂尤其對(duì)于內(nèi)點(diǎn)法。檢查模型規(guī)模如果可能?chē)L試從業(yè)務(wù)邏輯上簡(jiǎn)化模型減少不必要的變量或約束。6. 從練習(xí)題到實(shí)際項(xiàng)目思維轉(zhuǎn)變與經(jīng)驗(yàn)分享做了這么多練習(xí)題最終目的是為了應(yīng)用于實(shí)際。在實(shí)際項(xiàng)目中線性規(guī)劃的應(yīng)用遠(yuǎn)比課本習(xí)題復(fù)雜。首先數(shù)據(jù)準(zhǔn)備就是一大挑戰(zhàn)。練習(xí)題里的數(shù)據(jù)都是干凈的、給定的?，F(xiàn)實(shí)中原料價(jià)格、資源消耗系數(shù)、市場(chǎng)需求等數(shù)據(jù)都來(lái)自數(shù)據(jù)庫(kù)、報(bào)表或預(yù)測(cè)模型可能存在缺失、異常或不確定性。在構(gòu)建模型前必須進(jìn)行嚴(yán)格的數(shù)據(jù)清洗和驗(yàn)證。例如單位是否統(tǒng)一歷史消耗系數(shù)是否適用于未來(lái)我曾在做一個(gè)生產(chǎn)排程項(xiàng)目時(shí)因?yàn)檎`將“噸”和“公斤”的數(shù)據(jù)混用導(dǎo)致整個(gè)模型失效教訓(xùn)深刻。其次模型的維護(hù)和更新至關(guān)重要。一個(gè)線性規(guī)劃模型不是一勞永逸的。市場(chǎng)條件、生產(chǎn)工藝、公司政策都在變模型必須能方便地更新參數(shù)甚至調(diào)整結(jié)構(gòu)。這就要求我們?cè)诰幊虝r(shí)要有良好的習(xí)慣將數(shù)據(jù)f,A,b等與模型邏輯分離最好從外部文件如Excel、CSV讀取數(shù)據(jù)。這樣業(yè)務(wù)人員只需更新數(shù)據(jù)文件而無(wú)需觸碰代碼。% 良好的實(shí)踐從外部文件讀取數(shù)據(jù) data readmatrix(production_data.csv); f data(1, :); % 第一行是目標(biāo)系數(shù) A data(2:end-1, :); % 中間行是不等式約束矩陣 b data(end, :); % 最后一行是約束右側(cè)值 % ... 然后調(diào)用linprog最后結(jié)果的呈現(xiàn)和解釋需要技巧。你不能只給老板看一堆x1123.456這樣的數(shù)字。你需要將其轉(zhuǎn)化為業(yè)務(wù)語(yǔ)言“下個(gè)月建議將A生產(chǎn)線產(chǎn)能的60%用于生產(chǎn)產(chǎn)品X這樣可以預(yù)期提升毛利15%但需要提前儲(chǔ)備B原料20噸。” 同時(shí)一定要附上關(guān)鍵的前提假設(shè)和靈敏度分析結(jié)論“這個(gè)建議是基于產(chǎn)品Y價(jià)格穩(wěn)定在每單位100元的預(yù)測(cè)。如果其價(jià)格下跌超過(guò)5%方案需要重新評(píng)估?！本€性規(guī)劃是一個(gè)強(qiáng)大的工具但它只是一個(gè)工具。它的輸出質(zhì)量完全取決于輸入模型的質(zhì)量和對(duì)現(xiàn)實(shí)問(wèn)題的抽象能力。練習(xí)題教會(huì)我們?nèi)绾问褂眠@個(gè)工具而真正的藝術(shù)在于如何用它來(lái)刻畫(huà)和解決那些模糊、復(fù)雜、動(dòng)態(tài)的真實(shí)世界問(wèn)題。這需要不斷的練習(xí)、反思以及最重要的——與實(shí)際問(wèn)題深度碰撞。