基于多目標螢火蟲算法的供應鏈生產效率與穩(wěn)健性研究

舒 彤，高小琴，陳 收，汪壽陽，2

(1.湖南大學工商管理學院，長沙410082;2.中國科學院數學與系統(tǒng)科學研究院，北京100190)

基于多目標螢火蟲算法的供應鏈生產效率與穩(wěn)健性研究

舒 彤1，高小琴1，陳 收1，汪壽陽1，2

(1.湖南大學工商管理學院，長沙410082;2.中國科學院數學與系統(tǒng)科學研究院，北京100190)

文章針對供應鏈存在的不穩(wěn)定性，運用情景規(guī)劃法考慮不同種類的供應鏈中斷。在供應鏈網絡設計時考慮其穩(wěn)健性，對供應鏈中斷下的生產效率與穩(wěn)健性進行權衡，并采用多目標螢火蟲算法求解這個權衡問題，使結果具有連貫性。研究結果表明，目標函數最優(yōu)時的總成本反而比供應鏈最有效時的總成本低，產生了成本優(yōu)勢，企業(yè)可以根據生產效率與穩(wěn)健性權衡的近似帕累托前沿圖選擇合適的生產效率與穩(wěn)健性以促進企業(yè)的持續(xù)發(fā)展。

供應鏈中斷；多目標螢火蟲算法；生產效率；穩(wěn)健性

1 模型建立

對于供應鏈，以往文獻大多只考慮一種產品和一個制造中心，如Aviral(2011)[1]所涉及的供應鏈。這種設計不太符合實際，這里設計多種產品、多個制造中心的供應鏈。大多數企業(yè)擁有生產產品的制造中心，把產品運往各地的配送中心和對產品有需求的客戶區(qū)，這種三級供應鏈的模式也被諸多學者采用(Aviral，2011)[1]。因此這里設計了多種產品、多個制造中心的三級供應鏈，供應鏈包括位置固定的多產品制造中心a，潛在的配送中心b和位置固定的客戶區(qū)c，如圖1所示。

圖1 多種產品多個制造中心的三級供應鏈

鑒于Magdalini(2014)[2]的數學模型，本文在供應鏈設計時引入供應鏈的穩(wěn)健性，設計供應鏈節(jié)點和鏈接中斷的混合整數線性規(guī)劃模型如下。

1.1 參數設置

①索引

a:制造中心；b:配送中心；c:客戶區(qū)；N:產品種類；s:情景集合

②決策變量

制造中心是否供應配送中心和配送中心是否供應客戶區(qū)的二進制變量分別定義如下：

定義整數變量以描述多級供應鏈：

R(N，a):制造中心a生產產品N的數量；Q(N，a，b):從制造中心a到配送中心b產品N的數量；Q(N，b，c):從配送中心b到客戶區(qū)c產品N的數量

③需求參數

D(N，c):客戶區(qū)c對產品N的年需求量

④效率參數

ηe:供應鏈生產效率；ηr:供應鏈穩(wěn)健性

⑤成本參數

c(b，f):建立配送中心b時，每年所攤銷的固定成本；C(N，b，v):建立配送中心b時，每年攤銷的產品N的單位變動成本；C(N，a):制造中心a生產每單位N產品的生產成本；C(N，b，h):配送中心b所發(fā)生的產品N的單位持有成本；C(N，a，b):每單位產品N從制造中心a到配送中心b的單位運輸成本；C(N，b，c):每單位產品N從配送中心b到客戶區(qū)c的單位運輸成本；u(N):機會成本

⑥距離參數

G(a，b):從制造中心 a到配送中心 b的距離；H(b，c):從配送中心b到客戶區(qū)c的距離

⑦概率參數

p(s):情景s發(fā)生的概率

⑧中斷產品數量參數

q(N，s):情景s下N產品中斷的數量

1.2 約束條件

①網絡結構約束

制造中心、配送中心和客戶區(qū)之間所有相關網絡結構約束總結為：

式(1)表示如果制造中心a服務配送中心b，那么配送中心b至少供應某個客戶區(qū)。式(2)說明如果建立配送中心b，那么客戶區(qū)c可能由配送中心b供應，也有可能不由配送中心b供應。只有當制造中心a供應配送中心b時，制造中心a才能提供產品N給配送中心b，因此有約束條件(3)，其中k是一個適當的大數，可令k=10，0000，0000。同樣的約束條件可以應用到配送中心b和客戶區(qū)c之間，如式(4)所示。式(5)和式(6)是單一來源約束，是為了確保每個配送中心和每個客戶區(qū)分別只能由一個制造中心和一個配送中心供應。

②物料平衡約束

假設不存在庫存積累或損耗，物料平衡約束總結為：

式(7)說明制造中心a提供給所有配送中心的產品N的數量應當等于制造中心a生產的產品N的數量。同樣，制造中心提供給配送中心b的產品N的數量應當等于配送中心b提供給所有客戶區(qū)的產品N的數量，如式(8)所示。式(9)是為了確保每個客服區(qū)的需求得到滿足。

③非負約束

所有連續(xù)變量必須非負：

為有效縮減搜索空間，供應鏈生產效率和穩(wěn)健性必須非負：

1.3 目標函數

供應鏈的建立既要考慮生產效率也要兼顧穩(wěn)健性，因此目標函數設定為最大化供應鏈生產效率和最大化穩(wěn)健性兩個相互沖突的目標。供應鏈生產效率用運營成本來詮釋，供應鏈穩(wěn)健性用預期中斷成本來詮釋。

c(O)max:供應鏈最穩(wěn)健下的運營成本；c(O)min:供應鏈最有效時的運營成本；c(E)min:供應鏈最穩(wěn)健下的預期中斷成本；c(E)max:供應鏈最有效時的預期中斷成本。

目標函數中的運營成本包括基礎設施成本、生產成本、配送中心的物料持有成本和運輸成本。

式(17)為基礎設施成本，其之所以發(fā)生是由于配送中心的建立。這里假設配送中心的成本由固定成本和變動成本組成，其中固定成本是確定的，只要配送中心建立就存在，而變動成本則取決于配送中心產品N單位變動成本每年攤銷額乘以數量。假定制造中心的生產成本以單位成本的速率與產品的產量成正比，制造中心的總生產成本如式(18)所示。式(19)表示配送中心的物料持有成本，其與配送中心發(fā)生的總吞吐量成正比。式(20)和式(21)分別表示制造中心到配送中心的運輸成本和配送中心到客戶區(qū)的運輸成本，運輸成本是產品數量、距離和單位運輸成本的函數。一般情況下貨車是滿載貨物行駛的，因此運輸成本的規(guī)模經濟效應在此忽略不計。

綜上所述，運營成本的表達式為：

目標函數中的預期中斷成本用情景法定義。情景法是一個古老的概念，從最早記錄時間開始，人們就已經對未來很感興趣并且把情景法作為間接探索未來社會和制度的工具(Bradfield，2005)[3]。本文采用情景規(guī)劃法計算和分析供應鏈不同中斷情景下發(fā)生的預期中斷成本，情景分別設定為制造中心節(jié)點中斷或配送中心節(jié)點中斷和制造中心到配送中心的鏈接中斷或配送中心到客戶區(qū)的鏈接中斷。預期中斷成本用情景s發(fā)生的概率、情景s下N產品中斷的數量和產品N的單位邊際利潤的乘積來表示：

其中，p(s)為情景s發(fā)生的概率；u(N)是機會成本，即產品N的單位邊際利潤；q(N，s)為情景s下N產品中斷的數量：

所以，預期中斷成本的表達式為：

2 多目標螢火蟲算法

多目標螢火蟲算法可以同時考慮最大化供應鏈生產效率和最大化供應鏈穩(wěn)健性兩個相互沖突的目標，并且不同于單目標算法得出的離散點，可避免單目標算法陷入局部最優(yōu)。

2.1 螢火蟲算法

螢火蟲算法是由Yang(2008-2014)提出并不斷完善，是基于理想化的螢火蟲閃爍行為特征：(1)螢火蟲是雌雄皆宜的，所以螢火蟲會吸引其他螢火蟲而不管其性別是雌性還是雄性。(2)螢火蟲的吸引度與其亮度成正比，并且隨著距離的增加而減少。因此對于任何兩個閃爍的螢火蟲，不太亮的那個會向更亮的那個移動。如果對于一個特定的螢火蟲，沒有比其更亮的，那么這個螢火蟲會隨機移動；(3)螢火蟲的亮度取決于目標函數的值[4-10]。

對于最大化問題，螢火蟲的熒光亮度可以簡單地設定為與目標函數值成正比。在螢火蟲算法中，螢火蟲的相對熒光亮度和吸引度都影響螢火蟲的移動，這里需要對其進行定義。為簡單起見，可以假設一個螢火蟲的吸引度是由其熒光亮度決定的，而其熒光亮度又取決于目標函數。

然而，螢火蟲的吸引度還和螢火蟲之間的距離相關，隨著距離的不同而不同，因此定義螢火蟲的吸引度為：

其中，β0為r=0處的吸引度，即光源處的最大吸引度；γ表示光強吸收系數，用以模擬熒光在空中傳播逐漸衰減的特性，可設為常數；rij為位置分別處于xi和xj的任意兩個螢火蟲i和j之間的笛卡爾距離，。值得注意的是，以上定義的r不局限于歐幾里得距離，任何能有效解決優(yōu)化問題的測量都可以作為距離r。

螢火蟲i被亮度更高的螢火蟲 j吸引并向 j移動時的位置更新定義為：

式中，α為步長因子，通常是常數；rand為[0，1]上服從均勻分布的隨機因子；式中第二部分是吸引度所導致，第三部分為隨機擾動項，為避免過早陷進局部最優(yōu)解。

2.2 多目標螢火蟲算法

對于多目標優(yōu)化，一種方法是將所有的目標組合成一個單一目標，這樣就可以使用單目標優(yōu)化算法而不用太多修改。例如，Theofanis(2011)[11]以這種方式詳細研究了直接使用螢火蟲算法解決多目標優(yōu)化問題。另一種方法是拓展螢火蟲算法直接產生帕累托最優(yōu)前沿，通過擴展螢火蟲算法的基本思路，Yang(2013)[4]提出了多目標螢火蟲算法。

多目標螢火蟲算法的尋優(yōu)過程是：首先定義目標函數；然后初始化螢火蟲群，使其盡可能均勻分布在搜索空間，這可以通過使用抽樣技術來實現；規(guī)定可容忍誤差或最大迭代次數后，評估所有螢火蟲的亮度或目標函數值并對每一對螢火蟲進行比較；如果螢火蟲 j支配螢火蟲i，螢火蟲i就按公式(27)向螢火蟲j移動，移動之后如果i不滿足約束條件，那么會重新生成一個螢火蟲；如果一個螢火蟲不受其他任何螢火蟲支配，那么就把該螢火蟲放入帕累托前沿，生成一個隨機向量(之和等于1)，這樣就可以獲得一個最佳的組合解；接著非支配解集傳遞到下一次迭代；通過多次迭代，達到最大迭代次數后，一般可以得到近似帕累托前沿的n個非支配解集，從而實現尋優(yōu)。

為更有效地隨機移動，可通過最小化目標函數的加權和得到當前最優(yōu)解，此時：

從帕累托前沿的角度來看，如果一只螢火蟲沒有受到其他螢火蟲的支配，那么這只螢火蟲移動的位置更新為：

綜上所述，供應鏈中斷下生產效率與穩(wěn)健性權衡的多目標螢火蟲算法流程如下。

步驟1：定義目標函數，初始化螢火蟲群 xi(i=1，2，…，n)。

步驟3：通過非支配解更新拍累托前沿，記錄最優(yōu)解的個數，并把所有非支配解傳遞到下一次迭代，更新螢火蟲的亮度和位置。

步驟4：重復步驟2，直到達到最大的迭代次數，得到全部帕累托最優(yōu)解，找到當前最佳的近似帕累托前沿。

3 算例

為說明多目標螢火蟲算法對多種產品多個制造中心的三級供應鏈混合整數線性規(guī)劃模型的適用性，這里采用了一個位于中國的電腦制造公司。該公司生產多種產品，擁有多個制造中心和配送中心，符合多種產品、多個制造中心的三級供應鏈模型特征，因此用此算例來驗證多目標螢火蟲算法的適用性。該公司有兩個位于發(fā)達城市的制造中心，即北京和上海。公司選擇7個潛在的配送中心，分別位于華東、華南、華中、華北、西北、西南和東北地區(qū)交通最發(fā)達的省和直轄市，即北京、遼寧、上海、湖北、廣東、四川和陜西?？偣灿?3個客戶區(qū)，包括香港和澳門在內的33個大陸省級行政區(qū)。

3.1 問題闡述

假設兩個制造中心都能滿足所有需求，并且都生產兩種產品；一種為功能性產品：普通電腦，單位邊際利潤為200元；另一種為創(chuàng)新產品：新型電腦，單位邊際利潤為1300元。Fisher(1997)[12]根據產品的需求模式把產品分為兩類：即功能性產品和創(chuàng)新產品。功能性產品需求可預測，產品生命周期長，邊際貢獻低；創(chuàng)新產品的需求量不可預測，產品生命周期短，邊際貢獻高。例如，鹽、紙巾、牙刷等為功能性產品；新推出的汽車、時尚包等為創(chuàng)新產品。配送中心可以位于多達7個的潛在配送中心，最佳的配送中心位置由模型決定。制造中心、配送中心和客戶區(qū)之間的距離用百度地圖推薦路線計算。因為這項研究的重點是供應鏈節(jié)點中斷和鏈接中斷，因此這里假定客戶區(qū)的需求是確定的，并且與客戶區(qū)人口數成正比，省級行政區(qū)人口數根據《中國統(tǒng)計年鑒2013》得出[13]。

3.2 客戶區(qū)需求

33個客戶區(qū)對2種產品的年需求量如表1所示：普通電腦總的年需求量為135562臺，新型電腦總的年需求量為13558臺。

表1 客戶區(qū)對功能性產品和創(chuàng)新產品的需求

3.3 制造中心、配送中心和客戶區(qū)之間的距離

制造中心、配送中心和客戶區(qū)之間的距離用百度地圖推薦路線計算；從制造中心a到配送中心b的距離G(a，b)如表2所示，從配送中心b到客戶區(qū)c的距離H(b，c)如表3所示。

表2 制造中心與配送中心之間的距離(千米)

表3 配送中心與客戶區(qū)之間的距離(千米)

3.4 供應鏈中斷概率

供應鏈中斷可能由自然和人為因素造成，盡管自然災害發(fā)生的概率難以量化，但是歷史數據可用來預測自然災害發(fā)生的概率。Li(2013)[14]用歷史數據預測了未來自然災害發(fā)生的概率。這里假定在政治局勢穩(wěn)定的中國，大多數供應鏈中斷是由自然災害引起。為得出中斷的概率，這里采用《中國統(tǒng)計年鑒2013》中的分地區(qū)自然災害損失情況來計算省級行政區(qū)發(fā)生自然災害的相對概率[13]。制造中心和潛在配送中心中斷的概率如表4所示。

表4 制造中心以及配送中心中斷概率

3.5 相關成本

假設建立每個配送中心需要100000000元，且每個配送中心的使用年限為20年，因此建立配送中心時，固定成本每年所攤銷的費用為500000元。建立配送中心時，普通電腦單位變動成本每年攤銷費用為200元；新型電腦單位變動成本每年攤銷費用為800元。制造中心生產每臺普通電腦的生產成本都為2500元，制造中心生產每臺新型電腦的生產成本都為6000元。每個配送中心發(fā)生的普通電腦的單位持有成本均為50元，新型電腦的單位持有成本都為100元。每臺普通電腦和新型電腦從制造中心到配送中心的單位運輸成本都為40元，從配送中心到客戶區(qū)的單位運輸成本也為40元。對應算例在Matlab R2014b中使用多目標螢火蟲算法編程運行。

4 結果分析

供應鏈中斷可能是節(jié)點中斷，即配送中心中斷，也可能是鏈接中斷，即制造中心和配送中心之間的鏈接中斷。這里對這兩種情況都進行研究分析。

4.1 配送中心中斷

由于配送中心位置相距很遠，這里假定每個配送中心的中斷是相互獨立的并且多個配送中心的中斷可同時發(fā)生。每個配送中心只有兩種狀態(tài)：正常運作或中斷，并且假定配送中心中斷就會失去所有容量。配送中心中斷的概率取決于其所在的省級行政區(qū)，配送中心中斷的概率如表5所示。配送中心可能發(fā)生一個中斷，也可能多個中斷同時發(fā)生。這里最多考慮三個同時中斷，因為四個或四個以上配送中心同時中斷的概率很小。該模型具有1586個約束條件、739個變量和種情景。

在運行模型之前，c(O)的上下界和c(E)的上下界必須計算。c(O)min直接通過最小化c(O)得出，c(E)在這一點就是最大值。c(E)min可以直接通過最小化c(E)計算，此時選出的倉庫是最穩(wěn)健的，通過選出的倉庫求出此時的c(O)就是最大值。在配送中心中斷情況下，c(O)和c(E)的界限值如表5所示。

表5 配送中心中斷情況下c(O)和c(E)的界限值

模型使用多目標螢火蟲算法進行仿真，通過200次迭代后，配送中心中斷下最有效的供應鏈網絡如圖2所示。

圖2 配送中心中斷下最有效的供應鏈網絡

供應鏈最有效時的總成本(CO+CE)、最穩(wěn)健時的總成本和目標函數最優(yōu)時的最小總成本如表6所示。從表6可知，目標函數最優(yōu)時的總成本反而比供應鏈最有效時的總成本低，這是由于目標函數最優(yōu)時運營成本雖然增加了，但由此引起的預期中斷成本的減少額度比運營成本的增加額度要多。在市場競爭日益激烈的當今，成本降低會給企業(yè)帶來競爭優(yōu)勢。目標函數最優(yōu)總成本最小時的供應鏈網絡如圖3所示。對比圖2和圖3可知，北京配送中心轉移了一部分產品到湖北配送中心，由湖北配送中心配送，這是由于湖北配送中心的中斷概率比北京的中斷概率小所致。

圖3 配送中心中斷下目標函數最優(yōu)時的供應鏈網絡

使用多目標螢火蟲算法，經過200次迭代，獲得近似帕累托前沿，如圖4所示，橫軸表示供應鏈生產效率，縱軸表示供應鏈穩(wěn)健性。從圖4可知，供應鏈穩(wěn)健性隨供應鏈生產效率的增加而減小，企業(yè)可根據所處產業(yè)性質和自身情況選擇合適的供應鏈生產效率和穩(wěn)健性。比如，容易受自然災害影響的產業(yè)可適當提高供應鏈穩(wěn)健性，對自然災害不太敏感的產業(yè)可加大對生產效率的重視。

4.2 制造中心和配送中心之間的鏈接中斷

圖4 配送中心中斷下生產效率和穩(wěn)健性權衡的多目標帕累托前沿

鏈接是兩個節(jié)點之間的鏈接，因此假定鏈接中斷的概率為所連接節(jié)點中斷概率的平均值，制造中心和配送中心之間的鏈接中斷概率如表7所示。制造中心和配送中心之間的鏈接可能是一條中斷，也可能多條同時中斷，因為三條或三條以上鏈接中斷的概率很小。為避免問題過于復雜，這里最多考慮兩條同時中斷。該模型具有1586個約束條件、739個變量和種情景。

表7 制造中心和配送中心之間的鏈接中斷概率

在運行模型之前，c(O)的上下界和c(E)的上下界必須計算。c(O)min直接通過最小化c(O)得出，c(E)在這一點就是最大值。c(E)min可以直接通過最小化c(E)計算，此時的c(O)就是最大值。在制造中心和配送中心之間鏈接中斷這種情況下，c(O)和c(E)的界限值如表8所示。

表8 制造中心和配送中心鏈接中斷情況下c(O)和c(E)的界限值

模型使用多目標螢火蟲算法進行仿真，通過200次迭代后，制造中心和配送中心鏈接中斷下最有效的供應鏈網絡和目標函數最優(yōu)總成本最小時的供應鏈網絡分別如圖5和圖6所示。對比圖5和圖6可知，北京配送中心轉移了一部分產品到上海和陜西配送中心，由上海和陜西配送中心配送，這是由于上海和陜西配送中心的中斷概率比北京的中斷概率小所導致的。

圖5 制造中心和配送中心之間的鏈接中斷下最有效的供應鏈網絡

圖6 制造中心和配送中心中斷下目標函數最優(yōu)時的供應鏈網絡

使用多目標螢火蟲算法，經過200次迭代，獲得近似帕累托前沿，如圖7所示，橫軸表示供應鏈生產效率，縱軸表示供應鏈穩(wěn)健性。從圖4和圖7可知，無論節(jié)點中斷還是鏈接中斷，供應鏈穩(wěn)健性都隨供應鏈生產效率的增加而減小。根據圖7，企業(yè)可以對供應鏈生產效率和穩(wěn)健性進行權衡，考慮所處內部環(huán)境和外部環(huán)境選擇合適的供應鏈生產效率和穩(wěn)健性。

圖7 鏈接中斷下生產效率和穩(wěn)健性權衡的的多目標帕累托前沿

5 結論

從多種產品多個制造中心的三級供應鏈混合整數線性規(guī)劃模型的算例可知，無論是節(jié)點中斷還是鏈接中斷，供應鏈穩(wěn)健性隨供應鏈生產效率的增加而減小，而企業(yè)既希望提升生產效率又希望提升穩(wěn)健性，但這兩者此消彼長，所以必須對供應鏈中斷下的生產效率和穩(wěn)健性進行權衡，使目標函數最優(yōu)。由仿真結果可知，目標函數最優(yōu)時的總成本反而比供應鏈最有效時的總成本低，所以目標函數最優(yōu)時的近似帕累托前沿給企業(yè)提供了選擇合適供應鏈生產效率和穩(wěn)健性的依據?？偟膩碚f，這種方法可以給供應鏈規(guī)劃決策提供定量化工具。

[1]Shukla A，Lalit V A，Venkatasubramanian.Optimizing Efficiency-ro?bustness Trade-offs in Supply Chain Design under Uncertainty due to Disruptions[J].Physical Distribution and Logistics Management， 2011，41(6).

[2]Kalaitzidou M A，Longinidis P，Tsiakis P，et al.Optimal Design of Mul?tiechelon Supply Chain Networks With Generalized Production and Warehousing Nodes[J].Industrial and Engineering Chemistry Re?search，2014，53(33).

[3]Fielda R B，Wright G，Burt G，et al.The Origins and Evolution of Sce?nario Techniques in Long Range Business Planning[J].Futures，2005， 37(8).

[4]Yang X.Multiobjective Firefly Algorithm for Continuous Optimization [J].Engineering with Computers，2013，29(2).

[5]Yang X.Firefly Algorithms for Multimodal Optimization[C].Berlin: Springer，2009.

[6]Yang X.Firefly Algorithm，Lévy Flights and Global Optimization[C]. London:Springer，2010.

[7]Yang X，Deb.Eagle Strategy Using Lévy Walk and Firefly Algorithms for Stochastic Optimization[J].Studies in Computational Intelligence， 2010，284.

[8]Yang X.Firefly Algorithm，Stochastic Test Functions and Design Opti?misation[J].International Journal of Bio-inspired Compution，2012，2 (2).

[9]Yang X.Nature-Inspired Mateheuristic Algorithms:Success and New Challenges[J].Computer Engineering and Information Technology， 2012，1(1):1-3

[10]Yang X，Hosseini，Gandomi Amir Hossein.Firefly Algorithm for Solving Non-convex Economic Dispatch Problems With Valve Load?ing Effect[J].Applied Soft Computing，2012，12(3).

[11]Apostolopoulos T，Vlachos A.Application of the Fire fly Algorithm for Solving the Economic Emissions Load Dispatch Problem[J].Inter?national Journal of Combinatorics，Article ID 523806，2011.

[12]Marshall L.Fisher.What is the right supply chain for your product? [J].Harvard Business Review，1997，3(4).

[13]2013年中國統(tǒng)計年鑒[M].北京:中國統(tǒng)計出版社，2014.

[14]Li N，Liu X，Xie W，et al.The Return Period Analysis of Natural Di?sasters With Statistical Modeling of Bivariate Joint Probability Distri?bution[J].Risk Analysis，2013，33(1).

（責任編輯/易永生）

F406.2

A

1002-6487（2016）21-0049-06

國家自然科學基金資助項目(71172194；71221001；71390330；71390331)

舒 彤（1970—），男，江西波陽人，副教授，博士生導師，研究方向：供應鏈管理、商務智能。