□ 胡少龍 范婷睿

（西南交通大學 1.經濟管理學院，四川 成都 610031；2.服務科學與創(chuàng)新四川省重點實驗室，四川 成都 610031）

一、引 言

洪水、地震、颶風等自然災害易破壞基礎設施，造成大量災民無家可歸，需快速向災民提供應急物資。為了實現(xiàn)這一目標，需在災前選擇合適的區(qū)域儲備充足的物資，以便災后高效配送[1]。應急物資配置通常涉及設施選址和庫存兩個關鍵決策，其中設施選址決策包括物資儲備庫選址或臨時配送中心選址等；庫存決策指不同物資的儲備量。學者多應用混合整數(shù)規(guī)劃理論構建數(shù)學模型，優(yōu)化應急物資配置決策。劉波和李硯[2]、王蘇生等[3]構建了雙層規(guī)劃模型，優(yōu)化災后應急物資分配問題。葛春景等考慮多個受災點同時且多次需求的情況，構建了一個多重覆蓋選址的混合整數(shù)規(guī)劃模型，優(yōu)化物資配置決策[4]。于冬梅等從提高服務質量的視角，考慮設施中斷風險，建立了選址布局網(wǎng)絡的多目標混合整數(shù)規(guī)劃模型，優(yōu)化設施選擇決策[5]。

自然災害難以預測，受災范圍、發(fā)生時間等災情信息具有不確定性。為保障應急物資配置效率，需要考慮上述不確定因素。作為研究不確定環(huán)境下的最優(yōu)決策理論，隨機規(guī)劃是研究上述問題的有效工具。張慶和余淼[1]、 張夢玲等[6]、王海軍等[7]、Wang 等[8]構建隨機規(guī)劃模型，優(yōu)化災前物資配置和災后物資配送等決策。Sanci 和 Daskin[9]、Moreno 等[10]還考慮了道路受損的情況，構建隨機規(guī)劃模型，優(yōu)化應急物資配置和網(wǎng)絡恢復的聯(lián)合決策。Wang 等基于手機定位數(shù)據(jù)獲取災害信息構建情景，以優(yōu)化選址和配送等決策[11]。

面對大規(guī)模問題，優(yōu)化軟件和精確算法往往難以在可接受時間內求得最優(yōu)解，啟發(fā)式算法具有廣泛的應用場景。近年來，越來越多的學者將機器學習與優(yōu)化算法結合，以提高求解效率。李壯年等應用NSGA- Ⅱ遺傳算法求解多目標優(yōu)化問題，分別提出了基于特征工程和支持向量機等六種機器學習算法的參數(shù)優(yōu)化方法[12]。Fu 等應用統(tǒng)計機器學習理論降低隨機變量的維數(shù)和情景尺度[13]；Guevara 等提出了一種機器學習和分布魯棒優(yōu)化相結合的方法[14]；Ghasemi 等建立了基于機器學習的仿真元模型[15]；Zhang 等提出了一個基于極限學習機的集成學習模型和多目標規(guī)劃的優(yōu)化框架[16]。

樣本均值近似算法（Sample Average Approximation，SAA）是一種求解大規(guī)模隨機規(guī)劃問題的近似算法。Murali 等將禁忌搜索算法與SAA 結合[17]，Aydin 和Murat 提出了一種基于群體智能的SAA[18]，Li 和Zhang等在SAA 中引入了情景分解算法[19]，Bidhandi 和 Patrick 提出了一種基于加速采樣的方法改進SAA[20]，Jalali等提出了一種基于SAA 的遺傳算法[21]，Jiang 等提出了一種將SAA 與牛頓迭代相結合的方法[22]，Tao 等將SAA 與基于種群進化的人工藻類算法相結合[23]，以提高求解大規(guī)劃隨機規(guī)劃模型的效率。但是，以上研究難以保障生成最具代表性的樣本，以保證計算效率。如圖1(a) 所示，當一個樣本包含情景A 和C，另一個相比包含情景B 和C，則后者的代表性顯然較弱。這是由于情景B 和C 距離較近，所含有的信息比較相似。為解決該挑戰(zhàn)，Emelogu 等提出了一種基于機器學習的SAA，其中機器學習用于生成樣本[24]。作者利用聚類技術（如k-means）將相似情景進行聚類，然后隨機從每一簇選擇一個情景作為該簇的代表性情景，見圖1（b）和（c）。該算法的局限性在于，當每一簇包含的情景數(shù)量不同時，只選擇一個情景難以具有足夠的代表性。因此，本文的貢獻是提出整合分層隨機抽樣解決該挑戰(zhàn)。通過標準差決定每個簇的情景選取數(shù)量，可以改善樣本生成的合理性，以提高SAA 計算效率。

圖1 基于聚類的樣本生成過程

綜上所述，綜合考慮多種應急物資、不同的倉儲設施類別和需求的不確定性，以設施的選址、庫存和運輸成本，以及物資供應不足的懲罰成本最小為目標，設計決策變量分別表示倉儲設施選址、庫存和災后配送等決策，構建了基于情景的兩階段隨機規(guī)劃模型，優(yōu)化應急物資配置決策。此外，應用SAA 求解模型，并整合了一個機器學習框架以高效地生成樣本，進而提高SAA 的性能。所提出的機器學習框架使用k-means++ 對全部情景進行聚類，再應用分層隨機抽樣生成樣本。

其余部分內容如下。第二節(jié)，構建應急物資配置隨機規(guī)劃模型；第三節(jié)，設計基于機器學習的SAA 算法；第四節(jié)，設計數(shù)值實驗驗證算法的有效性；最后，對研究進行總結，并討論下一步研究方向。

二、應急物資配置隨機規(guī)劃模型

（一）問題描述

隨機規(guī)劃是一種處理優(yōu)化模型輸入值不確定性的技術之一，模型使用一組離散的情景表示不確定突發(fā)事件的影響范圍和影響程度。例如，一種情景表示雅安發(fā)生7級地震，影響了周圍市縣120 萬人；一種情景表示11 級臺風在寧波登陸，影響了周圍市縣50 萬人等等。兩階段隨機規(guī)劃是在不確定情況下，做出非預期的第一階段決策；第二階段決策是在第一階段決策和情景已知情況下進行[25]，即如何決定應急物資配置，使得應對各種災害情景的成本最小。具體地，構建情景集合描述不確定需求，以設施選址、采購、庫存、運輸和物資供應不足的懲罰成本最小的目標，優(yōu)化應急物資配置。應急物資配置決策包括：倉儲設施的選址和大小，每個設施中儲存的各種物資數(shù)量，以及應對不同情景時物資的配送、剩余和短缺數(shù)量。

（二）模型構建

1.模型假設

考慮飲用水、食物、醫(yī)療包三種應急物資；不同設施的儲備能力不同；不同物資的單位采購、運輸、持有和懲罰成本不同，且不隨情景變化。

2.符號說明

I 表示城市集合，由i,j索引。L表示設施類型集合，l∈L分別為大、中、小三種設施類型。A表示物資類別集合，a∈A分別為飲用水、食物、醫(yī)藥用品。S表示自然災害情景集合，s∈S。隨情景變化的參數(shù)包括情景發(fā)生概率sP和城市j對a類物資的需求Da,j,s，令未受災城市的需求量為0。儲備設施相關參數(shù)包括l類設施的容量Ul，l類設施的固定成本ClF。Hi,j表示城市i到j的距離，aV表示a類物資的單位體積，C aP表示a類物資的單位采購成本，C aT表示a類物資的單位運輸成本，CHa表示a類物資的單位持有成本，aG表示a類物資的單位懲罰成本。當情景s發(fā)生后，所儲備物資被迅速運往各受災城市，若庫存無法滿足受災城市的物資需求，則被視為物資供應不足。對此，設置懲罰成本aG，對未滿足物資進行懲罰。

第一階段決策變量包括：是否在i城市建造l型倉儲設施xi,l，取值分別為0 和1；i城市a類物資儲備量ya,i。第二階段決策變量包括：情景s下，i城市運送到j城市的a類物資的數(shù)量qa,i,j,s；情景s下，j城市的a類物資的短缺數(shù)量wa,j,s；情景s下，災害結束后i城市a類物資的剩余量za,i,s。

3.目標與約束函數(shù)

第一階段使設施固定成本和物資采購成本最小，即fc、pc最小。第二階段為當需求不確定時，運輸、倉儲和懲罰成本的期望最小。其中tc表示運輸成本，hc表示剩余物資的庫存成本，wc表示供應不足的懲罰成本。式（7）限制每個設施的采購總量不超過該設施的儲存能力，式（8）限制每個城市最多只能建設一種類別的倉儲設施，式（9）計算設施內應急物資的剩余數(shù)量，式（10）表示物資短缺量等于需求量減去已分配量。式（11）限定x只能取0或者1，式（12）為非負約束。

三、基于機器學習的樣本均值近似算法

樣本均值近似（SAA）是一種基于蒙特卡羅模擬的隨機離散優(yōu)化問題的求解方法。這種方法的基本思想是生成一個隨機樣本，然后用相應的樣本平均函數(shù)來近似期望值函數(shù)。對得到的樣本均值近似問題進行求解，該過程重復多次，直到估計目標值與下界的差值低于某一閾值時，停止迭代獲取滿意解[26]。SAA 中，更大的樣本意味著估計目標值更接近真實值。但是，隨著樣本的增大，求解的復雜度也隨之增長。因而，選擇合適的樣本量是SAA 的關鍵步驟。本節(jié)提出一種機器學習框架選取情景，以生成代表性樣本，提高SAA 求解效率。首先，應用k-means++ 聚類方法對災害情景進行聚類；然后，通過分層隨機抽樣方法獲得樣本。下面分別描述k-means++ 算法、分層隨機抽樣方法，以及改進后的SAA算法流程。

(一）k-means++算法

k-means算法也被稱為Lloyd 算法，由Stuart Lloyd 在1957 年提出。k-means方法是最流行的無監(jiān)督學習算法之一，它遵循一種簡單和相對有效的方法，將給定的數(shù)據(jù)集分類成一定數(shù)量的簇。首先，令K為中心集合，隨機選擇初始的個中心。然后，數(shù)據(jù)集中的每個點被分配到最接近它的中心的簇中。最后，計算每個簇的均值作為中心值，并根據(jù)新的中心值重復第二步，直到中心值不變則算法收斂。與k-means的隨機選擇初始中心相比，k-means++ 通過不同數(shù)據(jù)點之間的距離，迭代的選擇|K|個中心，以確保更快的收斂[27]。

k-means++算法步驟如下。

第一步：從需求數(shù)據(jù)集Da,j,s中隨機選取一個需求點作為聚類中心Ca,j,1。

第二步：通過輪盤賭方法確定下個中心點，分別計算需求點和已選中心的距離，并計算概率。

第五步：對每個簇gk，重新計算聚類中心

第六步：重復第4 和5 步，直到聚類中心位置不再變化。

(二）分層隨機抽樣算法

由于樣本是聚類后從各簇中抽取得到的，抽取方式對于樣本的選取也很重要。Emelogu 等應用簡單隨機抽樣，從每簇隨機抽取一個情景[24]，忽略了不同簇中情景數(shù)量差異較大時，選取一個情景難以具有代表性的問題。因此，提出了分層隨機抽樣算法應對該問題。首先，根據(jù)標準差計算各簇應抽取的情景數(shù)量。然后，通過隨機抽樣抽取情景。為防止隨機抽取的情景所對應的需求量過大或過小，通過限定其平均值，以確保樣本的合理性。例如，若抽取情景所對應的需求量都較大，模型偏好選取大型設施并儲備過多物資，會造成成本過大。

分層隨機抽樣算法步驟如下。

第一步：根據(jù)標準差計算各簇應抽取的情景數(shù)量[28]。

其中，σk為簇gk的標準差，表示樣本的大小，表示每個簇抽取的樣本數(shù)量；

第三步：若

重復第二步，其中γ∈{ 0,1}；

(三）算法對比

為對比分析Emelogu 等提出的算法[24]與本文提出的算法在求解大規(guī)模隨機規(guī)劃問題的效率，本節(jié)分別描述了兩個算法，應用胡少龍等[29]提出的SAA 框架。

本文提出的算法SKS：SAA，k-means++，以及分層隨機抽樣，如下。

第一步：應用k-means++算法；

第二步：應用分層隨機抽樣算法；

第三步：對任意m∈M，求解模型

第四步：使用全部情景集合S，且對任意m∈M，求解模型

Emelogu 等[24]提出的算法SKR: SAA，k-means++，以及簡單隨機抽樣，如下。

第一步：應用k-means++算法；

第二步：應用簡單隨機抽樣，每簇抽取一個情景；

第三步：同SKS 算法的第三步至第六步。

四、算例分析

（一）算例構造

構造颶風災害算例驗證算法的有效性。假設有一颶風多發(fā)區(qū)域，各城市受災情景受颶風登陸點和颶風等級影響，以下數(shù)據(jù)均來自已發(fā)表文獻。颶風通常用薩菲爾——辛普森級來劃分，Catg= {1 ,2,3,4,5}，其中5 級最嚴重。需求和概率兩個參數(shù)隨情景變化，其中，需求用受災人數(shù)衡量，其取決于受災城市的總人口、登陸城市和颶風等級。颶風的影響程度采用了Dalal 和üster 提出的方法表示[30]：

基于Rawls 和Turnquist 的研究[31]設立如下數(shù)據(jù)。考慮三種應急物資：水、食品和醫(yī)療包。假設水的單位是1000 加侖；食物為速食產品，以1000 餐為單位；醫(yī)療包為每人一個單位。表1 總結了應急物資的關鍵參數(shù)。假設每一種物資未滿足的懲罰成本為采購價格的10 倍，持有成本為購買價格的20%。大、中、小三種不同存儲容量的設施的固定成本和儲存容量，見表2。

表1 物資的單位采購價格、所占倉儲量和運輸成本

表2 不同類別設施的固定成本和儲存容量

（二）結果分析

表3 三種算法的計算時間和估計目標值

I，S M，N計算時間(秒) 目標值/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

為便于直觀比較，計算不同算例SKR和SKS的平均求解時間，見圖2。顯然，對于小規(guī)模算例，Gurobi具有一定的計算優(yōu)勢。而當數(shù)據(jù)量增大時（城市數(shù)量為40，且情景數(shù)量超過100），Gurobi求解時間陡然增大。SKR和SKS的計算時間增幅較小，計算效率明顯更高。此外，與SKR相比，當城市數(shù)量為20 時，SKS求解速度較快的優(yōu)勢并不明顯。但是，當城市數(shù)量增大到40 時，SKS的計算時間一直較少，表明其計算效率更高。當M不變，N值增大意味著情景數(shù)量變多，算法的求解時間主要取決于樣本問題的復雜性。例如，當I=40 、S=500、M=5 時，50 個情景的樣本問題顯然比25 個情景的更難求解。可通過整合更高效的算法求解樣本問題，以提高SKR和SKS的計算效率。當N不變，M值增大意味著樣本問題個數(shù)增多，求解時間則主要取決樣本問題的計算時間。例如，當I=40 、S=500、N=25 時，10 個樣本問題顯然要比5 個求解時間更長?？赏ㄟ^整合并行計算，以提高兩種算法的計算效率。

圖2 Gurobi、SKR 和SKS 平均計算時間的對比

以Gurobi求解結果為基準，計算Gap值并比較SKR和SKS的計算效率。如表3 所示，對于算例I=40 、S=500，Gurobi在3600 秒內無法求得最優(yōu)解，SKR和SKS可以在較短時間內找到更好的上界。因此，當問題規(guī)模增大到一定程度后，SKR和SKS在計算效率上都優(yōu)于Gurobi。但與SKR不同的是，無論M和N值為多少，SKS都能找到更好的上界。對于不同算例，與Gurobi求解結果相比，SKR的Gap在-0.57%至11.15%之間變化；而SKS僅在-0.56%至4.01%之間變化。為更直觀的對比不同算例下SKR和SKS的Gap值的變化，將不同算例下，兩個算法的結果求平均值，其結果見圖3。對于所有算例，顯然SKS所求得的上界都顯著優(yōu)于SKR。因此，通過對比分析計算時間和Gap值發(fā)現(xiàn)，SKS明顯優(yōu)于SKR。

圖3 SKR 和SKS 平均Gap 值的對比

五、結論

合理布局應急物資儲備庫并存儲物資可以在災害發(fā)生后及時做出響應，以縮短物資籌備及運輸?shù)臅r間。通過最大限度地提高應急物資供給能力，可以保障人民生命財產安全，降低突發(fā)事件所造成的損失。

考慮不確定需求，研究了應急物資配置的模型與算法，以優(yōu)化應急物資儲備庫選址與庫存決策。以設施選址、采購、庫存、運輸和供應不足的懲罰成本最小為目標，構建了應急物資配置的兩階段隨機規(guī)劃模型。同時，整合k-means++ 聚類算法和分層隨機抽樣，設計了一個機器學習框架生成樣本，以改進樣本均值近似方法，提高其在求解大規(guī)模隨機規(guī)劃問題的計算效率。最后，構建了災害情景對模型和算法進行了仿真驗證，通過與Gurobi和其他算法的對比分析表明，隨著算例規(guī)模增大，所提出的算法能在較短時間找到更好的上界。

盡管為求解大規(guī)模隨機規(guī)劃問題提供一條新的途徑，但研究仍存在不足，可從以下方面進一步完善。首先，僅考慮了單一的災害后果，未考慮到交通受阻、設施受損等的影響，需在后續(xù)的研究中加以討論。其次，如何有效確定簇的數(shù)量。k-means++算法容易理解，聚類效果好，尤其是在處理大數(shù)據(jù)集的時候，可以保證較好的伸縮性和高效率。但是，初始k 值需要人為設定，該數(shù)值選取常依靠經驗，可能帶來較大誤差。僅根據(jù)樣本數(shù)量設定k值，需要進一步改進。最后，未對比其他聚類算法的效果。下一步研究可以增加不同的聚類算法，通過實驗對比分析各聚類方法與SAA 結合后的優(yōu)劣。