不確定需求下容量限制工廠選址的改進蜂群算法

李睿雪，馬 良，劉 勇

（上海理工大學 管理學院，上海 200093）

隨著工業的不斷發展，工廠選址問題作為工業生產開展的第一環節，始終是工業中研究的重點。在以往的容量限制選址模型中通常將需求量作為確定值考慮，但是，在如今工業發展和市場環境不斷變化的情況下，需求點的需求不再是確定值。因此，考慮需求的不確定性對選址模型的影響顯得尤為重要。

近年來，國內外學者在選址問題中針對需求不確定以及容量限制作了大量的研究。在不確定需求方面，Bayram 等[1]在應急設施選址模型中考慮疏散需求的不確定性，建立了基于需求場景的兩階段隨機規劃疏散模型。張鑫等[2]通過區間灰數理論預測不確定需求，將差分進化和遺傳算法混合后對應急物資選址模型進行求解，但是，求解算法依賴大量的數據。成律等[3]在傳統配送中心選址中，利用隨機規劃理論和模糊理論構建了隨機客戶需求情況下配送中心選址的數學模型，并應用遺傳算法通過算例對模型進行驗證，但遺傳算法存在早熟等缺點。隨后國內外學者對容量限制因素進行研究，Jena等[4]發現容量限制對動態設施選址有極大影響。李建光等[5]通過離散場景表征避難需求的不確定性，考慮情形權重以及避難位置的規模，構建了一個魯棒隨機規劃模型來確定避難場所選址。商麗媛等[6]在多情景下需求不確定的應急物流配送中心選址模型中，將配送中心容量作為約束，對不確定需求用區間灰數表示，設計了免疫量子粒子群算法進行求解。徐小峰等[7]構建了考慮同時取送貨需求和帶模糊容量約束的在線設施動態選址模型。

然而，目前研究大多是考慮需求不確定或容量限制單一因素對工廠選址的影響，而實際中這兩種因素共同影響選址決策。基于此，本文建立不確定性需求下分級生產規模容量限制工廠選址模型。將備選工廠設置為按生產規模分級分段以滿足需求的不斷變化，并通過區間聯系數對不確定的需求進行預測。在算法上，之前的研究算法通常采用遺傳算法和粒子群算法等求解，存在早熟和收斂性差等缺點。由于本文提出的工廠選址模型考慮需求的不確定性和容量限制規模的變化，求解算法編碼和計算量相對之前研究更為復雜。而具有較好的適應性和并行性的蜂群算法更適用于本文模型。因此，采用人工蜂群算法求解以保證能夠達到全局最優解。

1 生產規模分階段容量限制工廠模型

1.1 工廠選址模型及參數

容量限制選址作為選址問題的擴展模型，其含義為在備選點滿足需求點需求并且不超過其本身的容量限制前提下，選擇總距離成本和建設成本最低的位置作為選址方案。針對需求不確定性這一特征，將生產規模設置為分階段規劃。根據備選位置不同需求量之和選擇合適的規模，更有利于節約成本以及后續發展。為此建立了改進的容量限制選址問題模型。

模型的選址目標是在備選位置中選擇能使距離成本最小和建設成本最小的最優位置。設I為備選工廠的位置集；J為需求點的位置集；K為時間集合；N為備選工廠不同生產階段規模的集合；ωj表示需求點j需求量歷年變化率的平均值；Rjk表示需求點j第k年的需求量變化率；dj為需求點j的需求量；tij表示從備選工廠位置i到需求點位置j之間的運輸成本；ein表示備選工廠位置i實際選擇的第n階段生產規模；gi表示備選工廠位置i的最大生產規模；hi表示備選工廠位置i在該地點的基礎建設成本；kin表示備選工廠位置i在第n階段生產能力下每單位的生產建設成本；f為選擇方案的總成本。

在工廠的生產規模大于需求點的需求之和的前提下，選擇最優的備選位置，模型如下：

模型中目標函數式（1）表示工廠選址的運輸成本和固定生產成本之和最小；式（2）表示需求點的需求量歷年變化率的平均值；約束式（3）表示備選工廠所服務的需求點的需求之和不超出該備選工廠的最大生產能力；式（4）表示確保需求點和開放的備選位置保持對應；式（5）和式（6）為決策變量約束，xij為工廠i對需求點j服務，yi為在備選工廠位置i選擇建設工廠，其中，開放取值為1，關閉取值為0。

1.2 基于聯系數的區間灰數的不確定需求預測

針對事物不確定性和確定性關系，本文提出聯系數這一概念，表示事物確定性測度和不確定性測度以及兩者之間關系。不確定需求的區間數通過聯系數表示，則可以更大程度地保留原來需求的不確定性[8]。

假設u=[A,B]是一個正區間數，可將其轉化為聯系數，表示為u=A+(B-A)i，根據聯系數定義，i∈[-1, 1]。轉化后區間數的范圍擴大了一倍，變為[A-(B+A),A+(B-A)]。但由于需求的實際情況，本文不需要負范圍，故取i∈[0, 1]。

灰色GM(1, 1)模型是常用的對含有不確定因素進行預測的模型，其中，灰色系統的基本元素為灰數，通常用“Gr”表示，指在某一個區間或者某個數集內取值的不確定數。區間灰數則是既有下界c又有上界d的灰數。記作Gr∈[c,d]，c≤d，上界與下界的差值為區間灰數的長度，記作l=d-c。在預測模型中需要將一些列區間灰數組成一個區間灰數序列，記作Gr{Gr1,Gr2, ···,Grn}[9]。

設區間灰數Gr∈[c,d]，則區間灰數序列為L{Gr}={Gr1,Gr2, ···,Grn}，可將區間灰數白化，其白化值為Gri=c+μd，μ=[0, 1]，μ為0～1的隨機數。同樣可以對區間灰數序列白化，其中，C{Gr}為白部，D{Gr}為區間序列白化后的灰部，如式（7）所示。

分別對區間灰色序列白部和灰部構建預測模型，將累加式還原到原始序列，得到式（8）和式（9）。

式中，β為最小二乘法的參數列中的參數。

結合聯系數和區間數之間的關系，將不確定需求的區間轉化為聯系數，再利用灰色預測模型將聯系數作為預測灰數，對需求進行預測。基于聯系數的區間灰數預測模型既適用于小樣本的數據，又能挖掘出樣本之間隱含的不確定性。

設區間灰數Gr∈[a,b]，U=A+Bi，A=a,B=b-a,i∈[0, 1]，則U=A+Bi為其對應的聯系數。U={U1,U2, ···,Un}，Uk=Ak+Bki，Ak=ak，Bk=bk-ak，其中，A和B分別為聯系數U的同部和異部，A={A1,A2, ···,An}和B={B1,B2, ···,Bn}為聯系數序列U的同步序列和異部序列。根據區間灰數預測模型分別對A和B構建預測模型，可得同部和異部的響應函數。

通過A和B值還原成區間灰數預測，即可得到需求量的模擬值。

2 改進的人工蜂群算法

2.1 人工蜂群算法

人工蜂群算法（artificial bee colony algorithm，ABC）是Karaboga等[10]在蜜蜂自組織采蜜的群體智能行為上提出的一種非數值優化計算方法，其中，蜂群實現采蜜的智能行為主要包含3個部分：蜜源、采蜜蜂、待工蜂，蜂群的采蜜過程主要包括：搜索蜜源、招募蜜蜂和放棄蜜源。

蜜源：蜜源表示存在的可能的解。通常為了方便求解，將收益度作為衡量蜜源優劣的標準；采蜜蜂：表示在蜜源中采蜜的蜜蜂，采蜜蜂通常攜帶蜜源信息，并且以一定比例成為引領蜂，引領蜂攜帶蜜源信息通過搖擺舞與其他蜜蜂分享這些信息；待工蜂：表示等待尋找蜜源采集的蜜蜂，他們可以分為兩種，偵察蜂和追隨蜂，偵查蜂搜索蜂巢附件的新蜜源，追隨蜂在蜂巢內等待,通過分享采蜜蜂的信息, 尋找蜜源。

2.2 人工蜂群算法的步驟

在ABC算法中，將蜜源的位置抽象為解空間中的點，即待求解問題的可行解。蜜源i的質量用收益度fit表示，np為蜜源的數量。ABC算法中引領蜂和待工蜂各占蜂群的一半，且同一時刻下同一蜜源中只有一只引領蜂采蜜[11]。

a. 設求解問題的維度為D，蜜源i的位置表示為X，其中，t為當前的迭代次數。蜜源的初始位置根據式（12）在搜索空間中隨機產生。

式中：xij代表第i個蜜源xi的第j維度值，i取值{1, 2, ···,sn}，j取值{1, 2, ···,d}；xmax,ij和xmin,ij分別為維度j的最大值和最小值。

初始化蜜源就是對每個蜜源的所在維度通過式（12）選取一個隨機值。

b. 引領蜂階段，在引領蜂階段，引領蜂用式（13）更新新蜜源。

式中：xkj代表領域蜜源，k取值{1,2, ···,sn}，且k不等于i；φ為取值在[-1, 1]之間的隨機數。

通過式（13）更新得到新蜜源，再借助貪婪算法的原則，比較更新前、后的蜜源適應度值，選擇最優者。

c. 追隨蜂階段，引領蜂搜索結束后，追隨蜂開始搜索，在這一階段，在舞蹈區接受引領蜂分享的蜜源信息，追隨蜂采用輪盤賭策略來選擇其所追隨的引領蜂。為了保證適應度更高的個體可以得到保留。同時對每個蜜源設置參數R，當蜜源更新后被保留，R為0；反之，則R加1。從而記錄一個蜜源沒有被更新的次數。

d. 偵察蜂階段。如果一個蜜源經過多次開采后，R值過高，并且超過了設定的閾值T，則需要放棄這個蜜源，同時進入偵察蜂階段。這一階段，偵察蜂利用式（12）隨機尋找新的蜜源，以代替被拋棄的蜜源。

e. 評估所有蜜蜂搜索得到食物源的適應度并進行排序，進行下一次迭代。

2.3 改進人工蜂群算法

基本的人工蜂群算法是為解決連續空間的問題而設計，解是由蜜蜂在蜜源附近搜索生成的，通常都是在一個連續的范圍空間。但是，在處理調度安排、選址方案等離散型問題時，蜜蜂的鄰域搜索從連續變為離散，這就使得算法原來的公式不再有效。為了使這一算法可以應用在離散問題上，本文基于人工蜂群算法的基本原理以及相關研究[12]，依照本文模型對公式進行改進，以適應離散化的選址模型，提出了改進的離散人工蜂群算法（new discrete artificial bee colony algorithm,NDABC），使得人工蜂群算法在保留原有并行性的優點下，對離散問題也有良好的適應性。同時也對其他算法的研究提供思路。

a. 個體編碼。建立蜜蜂個體的編碼形式是選址方案表示的關鍵，本文構建的工廠選址問題主要包括兩方面：每個備選點是否被選中和每個需求點和備選點的對應關系。本文將備選點按位置進行排序，設定標號（1, 2, ···,m），需求點同樣是按位置排序后形成一列n個元素的數列。

例如，有8個需求點，4個備選點，在4個備選點中選取2個位置作為工廠位置，解X為：

其中，X表示在4個備選位置中選取3和1號位置，3號位置為需求點1，3，4，6服務，1號位置為需求點2，5，7，8服務。

b. 初始化種群。作為整個算法計算的起點，優秀的起點可以更好地將蜜蜂分散在可行解中。為了更好地讓蜜蜂分布在可行范圍內，初始解采用隨機生成的辦法。備選點的選址數量和位置采用隨機生成方式，以擴大選址方案。在確定備選位置后，隨機在備選位置上分配需求點，并且將選址位置的方案作為蜜源，保證在初始生成的種群中個體分布在不同蜜源上。

為了保證蜜蜂群體的隨機搜索能力，在蜜蜂群每次迭代更新中將群體的一小部分待工蜂作為偵察蜂。通過數據實驗發現，根據適應度值進行排序后的45%個體作為追隨蜂；選取總體種群中5%的個體作為偵察蜂。這種比例分配可以達到較好的優化效果，保證進行隨機搜索，防止陷入局部最優。

c. 蜜蜂更新規則。

針對每一個個體都在自身的蜜源內進行搜索更新，本文參考了文獻[13]的更新策略，提出了一種多方式的領域更新策略。其中，包括3種領域更新策略：領域交換更新策略，記作U1，在待更新個體中隨機選取2個位置，并交換其值；領域插入更新策略，記作U2，在該更新個體中隨機選取2個位置，插入1個備選位置的序號，后面的元素依次后延；領域內倒轉更新策略，記作U3，待更新個體中隨機選取2個位置，將2個位置間的位置進行倒置更新。表1列出各策略的具體示例。

表1 蜜蜂領域更新策略Tab.1 Strategy of bee field renewal

d. 觀察蜂更新策略。

在人工蜂群算法中追隨蜂是在所追隨的采蜜蜂的蜜源附近進行搜索采蜜，所以，使用式（13）確定觀察蜂的蜜源Y=[ya,yb, ···]，其中，a,b表示備選位置的編號，Sf(i)表示在整個選址方案中選取備選位置i所產生的成本，同時保留其追隨的采蜜蜂的選址策略。

e. 追隨更新。在傳統蜂群算法中，蜜蜂是在不同蜜源間進行搜索，針對選址問題，首先將工廠位置按照地理位置進行排序，將選址方案中的備選點位置作為蜜源。因此，追隨蜂則會在所追隨的選址方案中隨機地在選取位置的前后進行變化，以模擬追隨蜂在蜜源附近搜索的行為。當追隨蜂表示的備選點確定后，則開始進行備選點與需求點間的對應。其中，更新原則按照上文的更新策略進行更新。

f. 偵察蜂更新策略。總體種群中選取總量5%的個體作為偵察蜂，進行隨機搜索，以廣大搜索范圍，防止陷入局部最優。

g. 局部搜索。在傳統人工蜂群算法中，為了防止單個蜜蜂多次在同一蜜源進行搜索，設置了最大搜索次數，在本算法中為了進一步擴大算法的搜索范圍，防止陷入局部最優。每個蜜蜂個體均在個體更新過程中更新R次，保留R次中適應度值最大的個體。

h. 選擇策略。采取類錦標賽，根據模型實際情況，與傳統錦標賽算法中選擇策略不同，本策略中在采蜜蜂個體中隨機選擇n個個體，找出這n個個體中的最大值，作為觀察蜂的追隨個體。這樣可以擴大隨機搜索方位，提高算法的搜索能力。

2.4 算法過程

a. 初始化蜜蜂種群，初始時刻所有蜜蜂都是偵察蜂，全局搜索蜜源并計算收益度；

b. 對步驟a中所有蜜蜂根據收益度的相對大小，蜜蜂轉化為采蜜蜂(收益度靠前)和觀察蜂，觀察蜂根據采蜜蜂的舞蹈來接收招募；

c. 每個采蜜蜂都會在原蜜源附近采蜜(局部搜索)；

d. 觀察蜂按照蜜源以一定概率選擇蜜源，并在附近局部搜索；

e. 采蜜蜂和觀察蜂都在其蜜源附近使用局部搜索策略，保留最優的個體，進入下一次更新；

f. 記錄當前所有蜜蜂找到的最優蜜源，執行步驟b。直至算法迭代結束。

3 仿真實驗與分析

為了驗證改進算法的有效性，以及模型的可行性，由于真實數據難以獲取，本文通過模擬文獻[14]的數據進行仿真實驗。依照文獻[14]的數據特征，通過Matlab在其基礎上模擬生成更大規模的數據，保證算例可以更好地模擬實際環境并能更好地驗證算法和模型。本文中備選點為20個，需求點為150，每個備選點可分為4階段的生產規模，每階段的單位生產成本為[16，14，12，10]（元/萬件），數據如表2所示。對每個需求點記錄近5年的需求，如表3所示。

表2 備選工廠位置參數Tab.2 Alternative plant location parameters

表3 需求點位置參數Tab.3 Location parameters of demand points

為了方便對比不同算法的迭代效果，本文加大了文中算法的迭代次數，減少群體的規模數量。保證得到更明細的算法優化結果曲線。因此，改進蜂群算法參數設置：初始種群個體數為40，迭代500次。同時通過遺傳算法（genetic algorithm, GA）、粒子群算法（particle swarm optimization, PSO）、海鷗算法（seagull optimization algorithm, SOA）和文獻[15]提出的改進螢火蟲算法（firefly algorithm, FA）作為對比組，改進蜂群算法參數設置：初始種群個體數為40，迭代500次；遺傳算法的初始種群為40，迭代500次，采用單點交叉和單點變異，其中交叉概率為0.5，變異概率為0.05；粒子群算法的初始種群為40，最大迭代次數500次，慣性因子為0.6，C1=C2=2，C1，C2為學習因子；螢火蟲算法的初始種群為40，迭代500次，其中，移動概率0.7，變異概率0.1。海鷗算法的初始種群為40，最大迭代次數為500次，控制因子Fk的值從2線性降低到0。分別對每個算法進行20次獨立計算后對比每種算法的平均值和標準差。每個需求點近5年的需求量使用式（10）和式（11）進行預測，將其代入算法中進行計算。選址工廠數為2。計算后算法迭代如圖1所示。

圖1 不同算法迭代對比圖Fig.1 Iterative comparison of different algorithms

根據圖1結果分析，本文算法相較于其他算法在迭代中不僅種群整體更新，而且加入了防止陷入局部最優的自身更新，由此提高了搜索效率，能夠快速穩定地得到最優解，迭代次數遠小于其他算法，收斂速度有了明顯的提高。另外，從表4中也可以看出，相對其他4種算法，本文提出的改進的人工蜂群算法，在20次獨立計算的平均值以及標準差明顯優于其他算法。在選址方案中，本文算法并結合其他算法得出的方案，1號位置及附近和12號位置及附近都是選址方案中較優的位置。本文提出的算法得到的選址方案剛好選擇1號和12號位置，并且在幾種算法對比中得到最優的總成本。根據備選工廠的數據可以看出，12號位置的生產規模比1號位置更大。本文的方案12號位置相對1號位置服務更多的需求點，可以更好地滿足需求點的需求。最后結合實際分析選址位置可以看出，本文算法所選的工廠位置更加分散，有利于實際情況下工廠的發展。

表4 不同算法計算結果對比Tab.4 Comparison of calculation results of different algorithms

在驗證了改進算法有效后，對模型本身進行深入分析。開始考慮改變容量限制工廠選址模型的選址數量對整個選址模型的影響，對于本文的工廠選址問題模型選址位置不易過多，因此，僅分別求解選址數量N=2，3，4這3種情況，每種情況獨立運行20次，計算20次的平均值、標準差以及迭代情況，運行結果如表5和圖2所示。

圖2 不同工廠選址數量的迭代曲線圖Fig. 2 Iterative curve of different numbers of plant sites

表5 不同選址數量下選址模型結果對比Tab.5 Comparison of location model results under different location numbers

通過分析不同選址數量的選址方案結果得出，當N=2時，即僅選擇2個備選位置作為工廠位置，得到的成本最小，所選擇的工廠的位置趨近于中心位置，其生產規模更大，服務需求點的數量更多，每個備選位置的后續擴建能力較差；而當N=3時，選擇3個備選位置產生的成本相對2個備選位置有所增加，工廠有選擇更小的生產規模的情況，可以兼顧后續需求增加后的擴建需求；最后N=4時，整個模型的選址策略所需的總成本最大，工廠的位置相比較其他策略更加分散，備選工廠僅需要選擇最小的生產規模就可以滿足需求。因此，如果選址方案的目標是追求低成本，那么可以選擇更少的選址數量，以達到更低的選址方案成本。但是，如果為了后續的發展以及更大程度地滿足需求的不斷增長，可以在選址方案中適當增加選址數量，以滿足后續發展中需求的增加以及服務尚未發掘的潛在需求點。

最后在選址位置上分析，3種不同選址數量的方案中均選擇到了2號位置附近、9號位置和12號位置附近。由此分析可得，這3個位置是后續選址決策中可以重點考慮分析的位置，為后續的選址方案提供參考。

4 結束語

分析了現有容量限制選址問題中缺少對需求不確定性考慮和容量限制約束多為固定值這兩部分。同時將這兩部分因素加入容量限制選址模型中，用聯系數和區間灰數將不確定的需求轉化為固定值以及將固定的容量限制約束設置為分階段，建立了全新的容量限制選址模型。并且采用了針對模型的特征而改進的人工蜂群算法對模型進行求解。通過與其他算法對比，驗證了算法在此類問題上更加高效。但本文的研究僅為單一產品，下一步可以拓展為多種產品下的選址問題。