□ 胡少龍 范婷睿
(西南交通大學(xué) 1.經(jīng)濟(jì)管理學(xué)院,四川 成都 610031;2.服務(wù)科學(xué)與創(chuàng)新四川省重點(diǎn)實(shí)驗(yàn)室,四川 成都 610031)
洪水、地震、颶風(fēng)等自然災(zāi)害易破壞基礎(chǔ)設(shè)施,造成大量災(zāi)民無家可歸,需快速向?yàn)?zāi)民提供應(yīng)急物資。為了實(shí)現(xiàn)這一目標(biāo),需在災(zāi)前選擇合適的區(qū)域儲(chǔ)備充足的物資,以便災(zāi)后高效配送[1]。應(yīng)急物資配置通常涉及設(shè)施選址和庫存兩個(gè)關(guān)鍵決策,其中設(shè)施選址決策包括物資儲(chǔ)備庫選址或臨時(shí)配送中心選址等;庫存決策指不同物資的儲(chǔ)備量。學(xué)者多應(yīng)用混合整數(shù)規(guī)劃理論構(gòu)建數(shù)學(xué)模型,優(yōu)化應(yīng)急物資配置決策。劉波和李硯[2]、王蘇生等[3]構(gòu)建了雙層規(guī)劃模型,優(yōu)化災(zāi)后應(yīng)急物資分配問題。葛春景等考慮多個(gè)受災(zāi)點(diǎn)同時(shí)且多次需求的情況,構(gòu)建了一個(gè)多重覆蓋選址的混合整數(shù)規(guī)劃模型,優(yōu)化物資配置決策[4]。于冬梅等從提高服務(wù)質(zhì)量的視角,考慮設(shè)施中斷風(fēng)險(xiǎn),建立了選址布局網(wǎng)絡(luò)的多目標(biāo)混合整數(shù)規(guī)劃模型,優(yōu)化設(shè)施選擇決策[5]。
自然災(zāi)害難以預(yù)測,受災(zāi)范圍、發(fā)生時(shí)間等災(zāi)情信息具有不確定性。為保障應(yīng)急物資配置效率,需要考慮上述不確定因素。作為研究不確定環(huán)境下的最優(yōu)決策理論,隨機(jī)規(guī)劃是研究上述問題的有效工具。張慶和余淼[1]、 張夢(mèng)玲等[6]、王海軍等[7]、Wang 等[8]構(gòu)建隨機(jī)規(guī)劃模型,優(yōu)化災(zāi)前物資配置和災(zāi)后物資配送等決策。Sanci 和 Daskin[9]、Moreno 等[10]還考慮了道路受損的情況,構(gòu)建隨機(jī)規(guī)劃模型,優(yōu)化應(yīng)急物資配置和網(wǎng)絡(luò)恢復(fù)的聯(lián)合決策。Wang 等基于手機(jī)定位數(shù)據(jù)獲取災(zāi)害信息構(gòu)建情景,以優(yōu)化選址和配送等決策[11]。
面對(duì)大規(guī)模問題,優(yōu)化軟件和精確算法往往難以在可接受時(shí)間內(nèi)求得最優(yōu)解,啟發(fā)式算法具有廣泛的應(yīng)用場景。近年來,越來越多的學(xué)者將機(jī)器學(xué)習(xí)與優(yōu)化算法結(jié)合,以提高求解效率。李壯年等應(yīng)用NSGA- Ⅱ遺傳算法求解多目標(biāo)優(yōu)化問題,分別提出了基于特征工程和支持向量機(jī)等六種機(jī)器學(xué)習(xí)算法的參數(shù)優(yōu)化方法[12]。Fu 等應(yīng)用統(tǒng)計(jì)機(jī)器學(xué)習(xí)理論降低隨機(jī)變量的維數(shù)和情景尺度[13];Guevara 等提出了一種機(jī)器學(xué)習(xí)和分布魯棒優(yōu)化相結(jié)合的方法[14];Ghasemi 等建立了基于機(jī)器學(xué)習(xí)的仿真元模型[15];Zhang 等提出了一個(gè)基于極限學(xué)習(xí)機(jī)的集成學(xué)習(xí)模型和多目標(biāo)規(guī)劃的優(yōu)化框架[16]。
樣本均值近似算法(Sample Average Approximation,SAA)是一種求解大規(guī)模隨機(jī)規(guī)劃問題的近似算法。Murali 等將禁忌搜索算法與SAA 結(jié)合[17],Aydin 和Murat 提出了一種基于群體智能的SAA[18],Li 和Zhang等在SAA 中引入了情景分解算法[19],Bidhandi 和 Patrick 提出了一種基于加速采樣的方法改進(jìn)SAA[20],Jalali等提出了一種基于SAA 的遺傳算法[21],Jiang 等提出了一種將SAA 與牛頓迭代相結(jié)合的方法[22],Tao 等將SAA 與基于種群進(jìn)化的人工藻類算法相結(jié)合[23],以提高求解大規(guī)劃隨機(jī)規(guī)劃模型的效率。但是,以上研究難以保障生成最具代表性的樣本,以保證計(jì)算效率。如圖1(a) 所示,當(dāng)一個(gè)樣本包含情景A 和C,另一個(gè)相比包含情景B 和C,則后者的代表性顯然較弱。這是由于情景B 和C 距離較近,所含有的信息比較相似。為解決該挑戰(zhàn),Emelogu 等提出了一種基于機(jī)器學(xué)習(xí)的SAA,其中機(jī)器學(xué)習(xí)用于生成樣本[24]。作者利用聚類技術(shù)(如k-means)將相似情景進(jìn)行聚類,然后隨機(jī)從每一簇選擇一個(gè)情景作為該簇的代表性情景,見圖1(b)和(c)。該算法的局限性在于,當(dāng)每一簇包含的情景數(shù)量不同時(shí),只選擇一個(gè)情景難以具有足夠的代表性。因此,本文的貢獻(xiàn)是提出整合分層隨機(jī)抽樣解決該挑戰(zhàn)。通過標(biāo)準(zhǔn)差決定每個(gè)簇的情景選取數(shù)量,可以改善樣本生成的合理性,以提高SAA 計(jì)算效率。
圖1 基于聚類的樣本生成過程
綜上所述,綜合考慮多種應(yīng)急物資、不同的倉儲(chǔ)設(shè)施類別和需求的不確定性,以設(shè)施的選址、庫存和運(yùn)輸成本,以及物資供應(yīng)不足的懲罰成本最小為目標(biāo),設(shè)計(jì)決策變量分別表示倉儲(chǔ)設(shè)施選址、庫存和災(zāi)后配送等決策,構(gòu)建了基于情景的兩階段隨機(jī)規(guī)劃模型,優(yōu)化應(yīng)急物資配置決策。此外,應(yīng)用SAA 求解模型,并整合了一個(gè)機(jī)器學(xué)習(xí)框架以高效地生成樣本,進(jìn)而提高SAA 的性能。所提出的機(jī)器學(xué)習(xí)框架使用k-means++ 對(duì)全部情景進(jìn)行聚類,再應(yīng)用分層隨機(jī)抽樣生成樣本。
其余部分內(nèi)容如下。第二節(jié),構(gòu)建應(yīng)急物資配置隨機(jī)規(guī)劃模型;第三節(jié),設(shè)計(jì)基于機(jī)器學(xué)習(xí)的SAA 算法;第四節(jié),設(shè)計(jì)數(shù)值實(shí)驗(yàn)驗(yàn)證算法的有效性;最后,對(duì)研究進(jìn)行總結(jié),并討論下一步研究方向。
隨機(jī)規(guī)劃是一種處理優(yōu)化模型輸入值不確定性的技術(shù)之一,模型使用一組離散的情景表示不確定突發(fā)事件的影響范圍和影響程度。例如,一種情景表示雅安發(fā)生7級(jí)地震,影響了周圍市縣120 萬人;一種情景表示11 級(jí)臺(tái)風(fēng)在寧波登陸,影響了周圍市縣50 萬人等等。兩階段隨機(jī)規(guī)劃是在不確定情況下,做出非預(yù)期的第一階段決策;第二階段決策是在第一階段決策和情景已知情況下進(jìn)行[25],即如何決定應(yīng)急物資配置,使得應(yīng)對(duì)各種災(zāi)害情景的成本最小。具體地,構(gòu)建情景集合描述不確定需求,以設(shè)施選址、采購、庫存、運(yùn)輸和物資供應(yīng)不足的懲罰成本最小的目標(biāo),優(yōu)化應(yīng)急物資配置。應(yīng)急物資配置決策包括:倉儲(chǔ)設(shè)施的選址和大小,每個(gè)設(shè)施中儲(chǔ)存的各種物資數(shù)量,以及應(yīng)對(duì)不同情景時(shí)物資的配送、剩余和短缺數(shù)量。
1.模型假設(shè)
考慮飲用水、食物、醫(yī)療包三種應(yīng)急物資;不同設(shè)施的儲(chǔ)備能力不同;不同物資的單位采購、運(yùn)輸、持有和懲罰成本不同,且不隨情景變化。
2.符號(hào)說明
I 表示城市集合,由i,j索引。L表示設(shè)施類型集合,l∈L分別為大、中、小三種設(shè)施類型。A表示物資類別集合,a∈A分別為飲用水、食物、醫(yī)藥用品。S表示自然災(zāi)害情景集合,s∈S。隨情景變化的參數(shù)包括情景發(fā)生概率sP和城市j對(duì)a類物資的需求Da,j,s,令未受災(zāi)城市的需求量為0。儲(chǔ)備設(shè)施相關(guān)參數(shù)包括l類設(shè)施的容量Ul,l類設(shè)施的固定成本ClF。Hi,j表示城市i到j(luò)的距離,aV表示a類物資的單位體積,C aP表示a類物資的單位采購成本,C aT表示a類物資的單位運(yùn)輸成本,CHa表示a類物資的單位持有成本,aG表示a類物資的單位懲罰成本。當(dāng)情景s發(fā)生后,所儲(chǔ)備物資被迅速運(yùn)往各受災(zāi)城市,若庫存無法滿足受災(zāi)城市的物資需求,則被視為物資供應(yīng)不足。對(duì)此,設(shè)置懲罰成本aG,對(duì)未滿足物資進(jìn)行懲罰。
第一階段決策變量包括:是否在i城市建造l型倉儲(chǔ)設(shè)施xi,l,取值分別為0 和1;i城市a類物資儲(chǔ)備量ya,i。第二階段決策變量包括:情景s下,i城市運(yùn)送到j(luò)城市的a類物資的數(shù)量qa,i,j,s;情景s下,j城市的a類物資的短缺數(shù)量wa,j,s;情景s下,災(zāi)害結(jié)束后i城市a類物資的剩余量za,i,s。
3.目標(biāo)與約束函數(shù)
第一階段使設(shè)施固定成本和物資采購成本最小,即fc、pc最小。第二階段為當(dāng)需求不確定時(shí),運(yùn)輸、倉儲(chǔ)和懲罰成本的期望最小。其中tc表示運(yùn)輸成本,hc表示剩余物資的庫存成本,wc表示供應(yīng)不足的懲罰成本。式(7)限制每個(gè)設(shè)施的采購總量不超過該設(shè)施的儲(chǔ)存能力,式(8)限制每個(gè)城市最多只能建設(shè)一種類別的倉儲(chǔ)設(shè)施,式(9)計(jì)算設(shè)施內(nèi)應(yīng)急物資的剩余數(shù)量,式(10)表示物資短缺量等于需求量減去已分配量。式(11)限定x只能取0或者1,式(12)為非負(fù)約束。
樣本均值近似(SAA)是一種基于蒙特卡羅模擬的隨機(jī)離散優(yōu)化問題的求解方法。這種方法的基本思想是生成一個(gè)隨機(jī)樣本,然后用相應(yīng)的樣本平均函數(shù)來近似期望值函數(shù)。對(duì)得到的樣本均值近似問題進(jìn)行求解,該過程重復(fù)多次,直到估計(jì)目標(biāo)值與下界的差值低于某一閾值時(shí),停止迭代獲取滿意解[26]。SAA 中,更大的樣本意味著估計(jì)目標(biāo)值更接近真實(shí)值。但是,隨著樣本的增大,求解的復(fù)雜度也隨之增長。因而,選擇合適的樣本量是SAA 的關(guān)鍵步驟。本節(jié)提出一種機(jī)器學(xué)習(xí)框架選取情景,以生成代表性樣本,提高SAA 求解效率。首先,應(yīng)用k-means++ 聚類方法對(duì)災(zāi)害情景進(jìn)行聚類;然后,通過分層隨機(jī)抽樣方法獲得樣本。下面分別描述k-means++ 算法、分層隨機(jī)抽樣方法,以及改進(jìn)后的SAA算法流程。
k-means算法也被稱為Lloyd 算法,由Stuart Lloyd 在1957 年提出。k-means方法是最流行的無監(jiān)督學(xué)習(xí)算法之一,它遵循一種簡單和相對(duì)有效的方法,將給定的數(shù)據(jù)集分類成一定數(shù)量的簇。首先,令K為中心集合,隨機(jī)選擇初始的個(gè)中心。然后,數(shù)據(jù)集中的每個(gè)點(diǎn)被分配到最接近它的中心的簇中。最后,計(jì)算每個(gè)簇的均值作為中心值,并根據(jù)新的中心值重復(fù)第二步,直到中心值不變則算法收斂。與k-means的隨機(jī)選擇初始中心相比,k-means++ 通過不同數(shù)據(jù)點(diǎn)之間的距離,迭代的選擇|K|個(gè)中心,以確保更快的收斂[27]。
k-means++算法步驟如下。
第一步:從需求數(shù)據(jù)集Da,j,s中隨機(jī)選取一個(gè)需求點(diǎn)作為聚類中心Ca,j,1。
第二步:通過輪盤賭方法確定下個(gè)中心點(diǎn),分別計(jì)算需求點(diǎn)和已選中心的距離,并計(jì)算概率。
第五步:對(duì)每個(gè)簇gk,重新計(jì)算聚類中心
第六步:重復(fù)第4 和5 步,直到聚類中心位置不再變化。
由于樣本是聚類后從各簇中抽取得到的,抽取方式對(duì)于樣本的選取也很重要。Emelogu 等應(yīng)用簡單隨機(jī)抽樣,從每簇隨機(jī)抽取一個(gè)情景[24],忽略了不同簇中情景數(shù)量差異較大時(shí),選取一個(gè)情景難以具有代表性的問題。因此,提出了分層隨機(jī)抽樣算法應(yīng)對(duì)該問題。首先,根據(jù)標(biāo)準(zhǔn)差計(jì)算各簇應(yīng)抽取的情景數(shù)量。然后,通過隨機(jī)抽樣抽取情景。為防止隨機(jī)抽取的情景所對(duì)應(yīng)的需求量過大或過小,通過限定其平均值,以確保樣本的合理性。例如,若抽取情景所對(duì)應(yīng)的需求量都較大,模型偏好選取大型設(shè)施并儲(chǔ)備過多物資,會(huì)造成成本過大。
分層隨機(jī)抽樣算法步驟如下。
第一步:根據(jù)標(biāo)準(zhǔn)差計(jì)算各簇應(yīng)抽取的情景數(shù)量[28]。
其中,σk為簇gk的標(biāo)準(zhǔn)差,表示樣本的大小,表示每個(gè)簇抽取的樣本數(shù)量;
第三步:若
重復(fù)第二步,其中γ∈{ 0,1};
為對(duì)比分析Emelogu 等提出的算法[24]與本文提出的算法在求解大規(guī)模隨機(jī)規(guī)劃問題的效率,本節(jié)分別描述了兩個(gè)算法,應(yīng)用胡少龍等[29]提出的SAA 框架。
本文提出的算法SKS:SAA,k-means++,以及分層隨機(jī)抽樣,如下。
第一步:應(yīng)用k-means++算法;
第二步:應(yīng)用分層隨機(jī)抽樣算法;
第三步:對(duì)任意m∈M,求解模型
第四步:使用全部情景集合S,且對(duì)任意m∈M,求解模型
Emelogu 等[24]提出的算法SKR: SAA,k-means++,以及簡單隨機(jī)抽樣,如下。
第一步:應(yīng)用k-means++算法;
第二步:應(yīng)用簡單隨機(jī)抽樣,每簇抽取一個(gè)情景;
第三步:同SKS 算法的第三步至第六步。
構(gòu)造颶風(fēng)災(zāi)害算例驗(yàn)證算法的有效性。假設(shè)有一颶風(fēng)多發(fā)區(qū)域,各城市受災(zāi)情景受颶風(fēng)登陸點(diǎn)和颶風(fēng)等級(jí)影響,以下數(shù)據(jù)均來自已發(fā)表文獻(xiàn)。颶風(fēng)通常用薩菲爾——辛普森級(jí)來劃分,Catg= {1 ,2,3,4,5},其中5 級(jí)最嚴(yán)重。需求和概率兩個(gè)參數(shù)隨情景變化,其中,需求用受災(zāi)人數(shù)衡量,其取決于受災(zāi)城市的總?cè)丝?、登陸城市和颶風(fēng)等級(jí)。颶風(fēng)的影響程度采用了Dalal 和üster 提出的方法表示[30]:
基于Rawls 和Turnquist 的研究[31]設(shè)立如下數(shù)據(jù)。考慮三種應(yīng)急物資:水、食品和醫(yī)療包。假設(shè)水的單位是1000 加侖;食物為速食產(chǎn)品,以1000 餐為單位;醫(yī)療包為每人一個(gè)單位。表1 總結(jié)了應(yīng)急物資的關(guān)鍵參數(shù)。假設(shè)每一種物資未滿足的懲罰成本為采購價(jià)格的10 倍,持有成本為購買價(jià)格的20%。大、中、小三種不同存儲(chǔ)容量的設(shè)施的固定成本和儲(chǔ)存容量,見表2。
表1 物資的單位采購價(jià)格、所占倉儲(chǔ)量和運(yùn)輸成本
表2 不同類別設(shè)施的固定成本和儲(chǔ)存容量
表3 三種算法的計(jì)算時(shí)間和估計(jì)目標(biāo)值
I,S M,N計(jì)算時(shí)間(秒) 目標(biāo)值/Gap Gurobi SKR SKS Gurobi(萬) SKR(%) SKS(%)3.61 4.01 10, 5 56 56 1.16 0.98 10, 10 62 58 3.63 0.31 5, 10 36 34 20, 100 19 6070 3.70 0.33 5, 20 58 62 11.15 0.35 10, 10 99 98 3.04 0.37 10, 20 110 102 4.01 0.19 5, 10 57 53 20, 200 40 6841 9.58 0.13 5, 50 143 122 9.20 0.00 5, 25 120 119 20, 500 135 7107 10, 25 224 240 8.33 0.08 10, 50 278 250 9.21 0.03 2.25 1.91 5, 10 152 110 3.90 0.17 10, 5 433 161 3.96 2.04 10, 10 310 190 3.87 0.17 5, 5 106 86 40, 100 92 12433 0.15 0.17 5, 20 633 449 1.04 0.15 10, 10 379 319 0.09 0.10 10, 20 734 800 0.57 0.18 5, 10 274 212 40, 200 3313 12092 4.67 -0.40 5, 50 1974 1068 -0.57 -0.56 5, 25 822 595 40, 500 3600 19025 10, 25 2154 1010 3.77 -0.44 10, 50 4808 3518 -0.57 -0.50
為便于直觀比較,計(jì)算不同算例SKR和SKS的平均求解時(shí)間,見圖2。顯然,對(duì)于小規(guī)模算例,Gurobi具有一定的計(jì)算優(yōu)勢(shì)。而當(dāng)數(shù)據(jù)量增大時(shí)(城市數(shù)量為40,且情景數(shù)量超過100),Gurobi求解時(shí)間陡然增大。SKR和SKS的計(jì)算時(shí)間增幅較小,計(jì)算效率明顯更高。此外,與SKR相比,當(dāng)城市數(shù)量為20 時(shí),SKS求解速度較快的優(yōu)勢(shì)并不明顯。但是,當(dāng)城市數(shù)量增大到40 時(shí),SKS的計(jì)算時(shí)間一直較少,表明其計(jì)算效率更高。當(dāng)M不變,N值增大意味著情景數(shù)量變多,算法的求解時(shí)間主要取決于樣本問題的復(fù)雜性。例如,當(dāng)I=40 、S=500、M=5 時(shí),50 個(gè)情景的樣本問題顯然比25 個(gè)情景的更難求解。可通過整合更高效的算法求解樣本問題,以提高SKR和SKS的計(jì)算效率。當(dāng)N不變,M值增大意味著樣本問題個(gè)數(shù)增多,求解時(shí)間則主要取決樣本問題的計(jì)算時(shí)間。例如,當(dāng)I=40 、S=500、N=25 時(shí),10 個(gè)樣本問題顯然要比5 個(gè)求解時(shí)間更長。可通過整合并行計(jì)算,以提高兩種算法的計(jì)算效率。
圖2 Gurobi、SKR 和SKS 平均計(jì)算時(shí)間的對(duì)比
以Gurobi求解結(jié)果為基準(zhǔn),計(jì)算Gap值并比較SKR和SKS的計(jì)算效率。如表3 所示,對(duì)于算例I=40 、S=500,Gurobi在3600 秒內(nèi)無法求得最優(yōu)解,SKR和SKS可以在較短時(shí)間內(nèi)找到更好的上界。因此,當(dāng)問題規(guī)模增大到一定程度后,SKR和SKS在計(jì)算效率上都優(yōu)于Gurobi。但與SKR不同的是,無論M和N值為多少,SKS都能找到更好的上界。對(duì)于不同算例,與Gurobi求解結(jié)果相比,SKR的Gap在-0.57%至11.15%之間變化;而SKS僅在-0.56%至4.01%之間變化。為更直觀的對(duì)比不同算例下SKR和SKS的Gap值的變化,將不同算例下,兩個(gè)算法的結(jié)果求平均值,其結(jié)果見圖3。對(duì)于所有算例,顯然SKS所求得的上界都顯著優(yōu)于SKR。因此,通過對(duì)比分析計(jì)算時(shí)間和Gap值發(fā)現(xiàn),SKS明顯優(yōu)于SKR。
圖3 SKR 和SKS 平均Gap 值的對(duì)比
合理布局應(yīng)急物資儲(chǔ)備庫并存儲(chǔ)物資可以在災(zāi)害發(fā)生后及時(shí)做出響應(yīng),以縮短物資籌備及運(yùn)輸?shù)臅r(shí)間。通過最大限度地提高應(yīng)急物資供給能力,可以保障人民生命財(cái)產(chǎn)安全,降低突發(fā)事件所造成的損失。
考慮不確定需求,研究了應(yīng)急物資配置的模型與算法,以優(yōu)化應(yīng)急物資儲(chǔ)備庫選址與庫存決策。以設(shè)施選址、采購、庫存、運(yùn)輸和供應(yīng)不足的懲罰成本最小為目標(biāo),構(gòu)建了應(yīng)急物資配置的兩階段隨機(jī)規(guī)劃模型。同時(shí),整合k-means++ 聚類算法和分層隨機(jī)抽樣,設(shè)計(jì)了一個(gè)機(jī)器學(xué)習(xí)框架生成樣本,以改進(jìn)樣本均值近似方法,提高其在求解大規(guī)模隨機(jī)規(guī)劃問題的計(jì)算效率。最后,構(gòu)建了災(zāi)害情景對(duì)模型和算法進(jìn)行了仿真驗(yàn)證,通過與Gurobi和其他算法的對(duì)比分析表明,隨著算例規(guī)模增大,所提出的算法能在較短時(shí)間找到更好的上界。
盡管為求解大規(guī)模隨機(jī)規(guī)劃問題提供一條新的途徑,但研究仍存在不足,可從以下方面進(jìn)一步完善。首先,僅考慮了單一的災(zāi)害后果,未考慮到交通受阻、設(shè)施受損等的影響,需在后續(xù)的研究中加以討論。其次,如何有效確定簇的數(shù)量。k-means++算法容易理解,聚類效果好,尤其是在處理大數(shù)據(jù)集的時(shí)候,可以保證較好的伸縮性和高效率。但是,初始k 值需要人為設(shè)定,該數(shù)值選取常依靠經(jīng)驗(yàn),可能帶來較大誤差。僅根據(jù)樣本數(shù)量設(shè)定k值,需要進(jìn)一步改進(jìn)。最后,未對(duì)比其他聚類算法的效果。下一步研究可以增加不同的聚類算法,通過實(shí)驗(yàn)對(duì)比分析各聚類方法與SAA 結(jié)合后的優(yōu)劣。