基于碳稅政策的整車物流網(wǎng)絡(luò)優(yōu)化模型及算法研究

（中南大學(xué)交通運(yùn)輸工程學(xué)院，長沙 410075）

引 言

在 《車輛運(yùn)輸車治理工作方案》出臺后，整車物流領(lǐng)域發(fā)生重大變革，自2018年7月1日起，不合規(guī)車輛運(yùn)輸車被全面禁止通行，整車公路運(yùn)輸成本急劇上升，企業(yè)開始采用圖1所示的整車“庫前移”物流網(wǎng)絡(luò)代替原有的公路直達(dá)運(yùn)輸網(wǎng)絡(luò)。該網(wǎng)絡(luò)涉及配送中心選址、庫存及運(yùn)輸?shù)拳h(huán)節(jié)，需求量的變化會導(dǎo)致這3個環(huán)節(jié)的成本變化，且由于整車需求量具有明顯的季節(jié)性規(guī)律，在對整車物流網(wǎng)絡(luò)進(jìn)行優(yōu)化時需考慮不同時期的不同需求量要求，使得決策方案具有較強(qiáng)的魯棒性。此外，由于環(huán)境問題日益嚴(yán)重，而物流配送作為碳排放的主要來源之一，考慮碳稅政策對整車物流網(wǎng)絡(luò)優(yōu)化方案的影響具有非常重要的現(xiàn)實(shí)意義。

圖1 整車物流網(wǎng)絡(luò)

針對低碳視角下的物流網(wǎng)絡(luò)優(yōu)化問題，Yang等［1］提出了碳稅約束型城市物流配送網(wǎng)絡(luò)規(guī)劃模型；Jin等［2］對不同碳減排政策下大型零售商的供應(yīng)鏈網(wǎng)絡(luò)優(yōu)化及運(yùn)輸方式選擇進(jìn)行了研究；成耀榮和譚維［3］構(gòu)建了基于低碳經(jīng)濟(jì)視角的多式聯(lián)運(yùn)路徑優(yōu)化模型；Saxena等［4］研究了碳稅和減排獎勵政策對印度輪胎再制造物流網(wǎng)絡(luò)設(shè)計的影響。針對整車物流網(wǎng)絡(luò)優(yōu)化問題，王沛等［5］研究了多層級、多周期的整車鐵路物流基地選址；秦璐等［6］利用正交試驗(yàn)方法對整車物流網(wǎng)絡(luò)優(yōu)化進(jìn)行了研究；鄭燕等［7］對長江流域整車公水聯(lián)運(yùn)網(wǎng)絡(luò)進(jìn)行了優(yōu)化。在物流網(wǎng)絡(luò)魯棒性研究方面，Xi等［8］研究了非等區(qū)域動態(tài)設(shè)施布局，提出了考慮貨物裝卸點(diǎn)位置的混合魯棒優(yōu)化模型；周愉峰等［9］針對國家血液戰(zhàn)略儲備庫選址，構(gòu)建了考慮多情景、多血型的隨機(jī)p-魯棒選址-庫存優(yōu)化模型；孫華麗等［10］考慮應(yīng)急需求的隨機(jī)性，建立了多目標(biāo)選址-路徑魯棒優(yōu)化模型。

已有文獻(xiàn)大多是單獨(dú)考慮物流網(wǎng)絡(luò)魯棒性或碳排放而建立的優(yōu)化模型，鮮有文獻(xiàn)同時考慮這兩者的共同影響并針對整車物流行業(yè)背景進(jìn)行研究。因此，本文同時考慮碳稅政策和多種需求量情景，建立考慮碳稅政策的整車物流網(wǎng)絡(luò)魯棒優(yōu)化模型，設(shè)計遺傳退火算法進(jìn)行求解，為整車物流網(wǎng)絡(luò)優(yōu)化問題提供更全面的解決思路。

1 問題描述與建模

1.1 問題描述

由主機(jī)廠、配送中心和4s店構(gòu)成的整車物流網(wǎng)絡(luò)，如圖1所示。由于不同時期整車的需求量具有較大差異，研究采用 “情景”處理方式，即設(shè)置若干個典型的需求情景，每種情景對應(yīng)不同的需求量。要解決的問題是：在碳稅政策下，確定配送中心的數(shù)量和選址位置、配送中心與4s店的服務(wù)關(guān)系、配送中心的庫存控制參數(shù)以及主機(jī)廠至配送中心的運(yùn)輸方式，使得該方案在滿足魯棒性約束的要求下，在所有需求量情景下的總成本期望值最小化。

1.2 模型假設(shè)

（1）配送中心容量能夠滿足需求；

（2）1個主機(jī)廠和若干4s店的位置已確定，整車從主機(jī)廠至4s店經(jīng)過且只經(jīng)過1個配送中心，不考慮從主機(jī)廠至4s店的直達(dá)運(yùn)輸；

（3）各情景發(fā)生的概率可以主觀設(shè)定；

（4）各個4s店的需求相互獨(dú)立，且單位時間 （天）內(nèi)的需求量均服從正態(tài)分布；

（5）主機(jī)廠至配送中心的干線運(yùn)輸可選擇公路、鐵路和水路運(yùn)輸中的一種，但從配送中心至4s店只采用公路運(yùn)輸進(jìn)行配送。

（6）配送中心采用連續(xù)檢查的（Q，R）庫存控制策略，而4s店只保有少量的庫存，忽略4s店的庫存成本。

1.3 模型參數(shù)及含義

I配送中心候選點(diǎn)集合，i∈I；

J4s店集合，j∈J；

O整車主機(jī)廠所在地點(diǎn)；

G需求情景集合，g∈G；

K運(yùn)輸方式集合 （k∈K，k＝1，2，3分別表示公路運(yùn)輸、鐵路運(yùn)輸和水路運(yùn)輸）；

ck運(yùn)輸方式k的單位運(yùn)輸成本 （元／（km·輛））；

ek運(yùn)輸方式k的碳排放因子（kg／（km·輛））；

doi主機(jī)廠o與配送中心i的距離；

dij配送中心i與4s店j的距離 （km）；

α庫存服務(wù)水平；

Zα庫存服務(wù)水平為α?xí)r的安全庫存系數(shù)；

si配送中心i的安全庫存量 （輛）；

hi配送中心i的儲存成本 （元／（輛·天））；

Ri配送中心i的訂貨點(diǎn) （輛）；

Ai配送中心i的平均庫存 （輛）；

TPi配送中心i連續(xù)兩次訂購的間隔時間（天）；

Q配送中心i的訂購批量 （輛）；

fi配送中心i的固定設(shè)施日折舊成本 （元／天）；

從主機(jī)廠至配送中心i采用運(yùn)輸方式k的提前期 （天）；

OC配送中心補(bǔ)貨的固定成本 （元／次）；

φ碳稅 （元／kg）；

qg情景g發(fā)生的概率；

p遺憾值限定系數(shù)。

決策變量：

Xi在節(jié)點(diǎn)i建立配送中心時為1，否則為0；Yijj由配送中心i服務(wù)時為1，否則為0；Uik從主機(jī)廠o至配送中心i采用運(yùn)輸方式k時為1，否則為0。

1.4 物流網(wǎng)絡(luò)總成本分析

（1）庫存成本分析

整車屬于高價值商品，需要經(jīng)常性的進(jìn)行庫存盤點(diǎn)，因此配送中心采用連續(xù)檢查的（Q，R）庫存控制策略。根據(jù)模型假設(shè) （4），由數(shù)理統(tǒng)計知識可知配送中心的需求服從正態(tài)分布，其中：

根據(jù)庫存控制相關(guān)理論，配送中心的日常庫存控制參數(shù)如下：

①提前期：

②安全庫存：

③訂貨點(diǎn)：

④平均庫存：

⑤連續(xù)兩次訂購的間隔時間：

將式 （8）對Qi進(jìn)行求導(dǎo)，并令其等于0，可得到配送中心i的最佳訂購批量， 如式 （9）所示：

將式 （9）代入式 （7）中即可確定TPi的值。

（2）運(yùn)輸成本分析

運(yùn)輸成本包括主機(jī)廠至配送中心的長途干線運(yùn)輸成本、鐵路及水路運(yùn)輸兩端的接駁換裝成本、以及從配送中心至4s店的配送成本，則單位時間內(nèi)的平均運(yùn)輸成本如式 （10）所示：

（3）碳排放成本分析

物流網(wǎng)絡(luò)中的總碳排放量包括干線運(yùn)輸、中轉(zhuǎn)以及配送產(chǎn)生的碳排放成本，則單位時間內(nèi)的平均碳排放成本如式 （11）所示：

（4）固定設(shè)施折舊成本分析

配送中心i的固定設(shè)施日折舊成本如式 （12）所示：

綜上所述，在情景g下系統(tǒng)總成本Zg可表示為各項成本之和，如式 （13）所示：

1.5 模型建立

（1）基于特定情景的確定性模型

對于某種特定情景g，其需求量參數(shù)是確定的，綜合以上的分析，可建立如下確定性優(yōu)化模型A：

式 （14）為最小化系統(tǒng)總成本，式 （15）表示每個4s店只能分配給1個配送中心；式 （16）表示干線運(yùn)輸只能選擇1種運(yùn)輸方式；式 （17）確保4s店不會分配給未選中的配送中心，式 （18）表示至少有1個備選配送中心被選中，式 （19）～（21） 為0～1變量約束。

（2）魯棒優(yōu)化模型

由模型A雖然可求得每個情景所對應(yīng)的優(yōu)化方案及目標(biāo)函數(shù)值Z?g，但由于不同情景下的需求量有較大差異，該情景下的優(yōu)化方案在其他情景中可能表現(xiàn)會很差。因此，在模型A的基礎(chǔ)上建立魯棒優(yōu)化模型B。模型B力求得到一個在所有情景下的目標(biāo)函數(shù)值與該情景下的最優(yōu)目標(biāo)函數(shù)值的差值在可接受范圍內(nèi)的滿意解，其定義如下所示［11］： 設(shè)為情景g下的最優(yōu)目標(biāo)函數(shù)值，r為滿足所有情景g∈G下的可行解，為r在情景g下的目標(biāo)函數(shù)值，如果解r滿足魯棒約束條件，則稱r為該問題的魯棒解，其中參數(shù)p被稱為遺憾值限定系數(shù)，反映了模型的魯棒性水平，p值越小，則模型的魯棒性越強(qiáng)。根據(jù)定義可知，滿足約束條件的魯棒解r有多個，其中必有一個最佳的解r?使得總成本最小。模型B如下所示：

其他約束條件同式 （15）～（21）， 式 （22） 為最小化魯棒解在各情景下成本的期望值，式 （23）表示各情景發(fā)生的概率之和為1，式 （24）表示魯棒約束條件。此外，當(dāng)不考慮約束條件 （24）時，模型B即等同于隨機(jī)規(guī)劃模型C。

2 模型求解算法

選址—庫存問題已被證明為NP難問題［12］。而本文所研究的整車物流網(wǎng)絡(luò)優(yōu)化涉及選址、庫存、運(yùn)輸方式選擇等，相比選址-庫存問題而言更為復(fù)雜，難以用精確算法進(jìn)行求解。在求解該類問題的啟發(fā)式算法中，遺傳算法的全局尋優(yōu)能力強(qiáng)，收斂速度快［13］，但其局部搜索能力較差，易產(chǎn)生 “早熟”現(xiàn)象［14］。而模擬退火算法具有較強(qiáng)的局部搜索能力，因此，本文結(jié)合這兩種算法優(yōu)點(diǎn)，設(shè)計遺傳退火算法對該問題進(jìn)行求解，即在遺傳操作陷入局部收斂后 （算法所得的歷史最優(yōu)解在一定代數(shù)內(nèi)不發(fā)生變化時可視為發(fā)生局部收斂），再引入模擬退火操作使算法跳出局部最優(yōu)解。先給出求解模型A的具體步驟，如下所示：

Step1：編碼。采用實(shí)數(shù)編碼，每條染色體的編碼分為兩段，第1段編碼包含n位基因，對應(yīng)4s店與配送中心的分配關(guān)系，取值范圍為［1，m］之間的整數(shù)；第2段編碼包含m位基因，對應(yīng)主機(jī)廠至配送中心的運(yùn)輸方式選擇變量，取值范圍為［1，3］之間的整數(shù)。

Step2：隨機(jī)生成規(guī)模為popsize的初始種群1，設(shè)置算法最大迭代次數(shù)為maxgen，當(dāng)前迭代次數(shù)為gen，令gen＝0。以目標(biāo)函數(shù)作為適應(yīng)度函數(shù)，記錄種群1中最優(yōu)個體的適應(yīng)度值作為歷史最優(yōu)解。

Step3：選擇。結(jié)合 “輪盤賭”與 “精英保留”策略對種群進(jìn)行篩選。

Step4：交叉。采用單點(diǎn)交叉方式，如圖2所示。

Step5：變異。采用交換變異方式，由于染色體的兩段編碼的取值范圍不同，為了避免變異后產(chǎn)生非法個體，只有當(dāng)隨機(jī)選擇的兩個變異點(diǎn)位于同一段編碼時才交換兩點(diǎn)的值。如圖2所示：

圖2 交叉、變異過程

Step6：記交叉變異后的種群為種群2。用種群1中的最優(yōu)個體替換種群2的最差個體，產(chǎn)生新種群3，計算適應(yīng)度，更新歷史最優(yōu)解及其保持不變的代數(shù)。

Step7：判斷是否滿足模擬退火操作執(zhí)行條件。當(dāng)歷史最優(yōu)解保持不變的代數(shù)超過100代時為模擬退火操作執(zhí)行的必要條件。若滿足該條件，則從種群3中隨機(jī)選擇比例為?的個體構(gòu)成模擬退火操作的初始解集，執(zhí)行Step8，否則，跳轉(zhuǎn)至Step13。

Step8：初始化模擬退火操作。令當(dāng)前解為w0，初溫為T0，溫度下降系數(shù)為low，每一溫度下的最大迭代次數(shù)為maxnum，令當(dāng)前溫度T＝T0，當(dāng)前溫度下已迭代次數(shù)num＝0。

Step9：隨機(jī)擾動生成新解w1。新解產(chǎn)生機(jī)制為隨機(jī)選取w0中染色體π的某個位置γ的基因作為擾動點(diǎn)，當(dāng)γ位于第1段編碼時，擾動范圍為［1，m］之間的整數(shù)，如式 （25）；當(dāng)γ位于第2段編碼時，擾動范圍為［1，3］之間的整數(shù)，如式（26）。其中ceil為MATLAB中朝正無窮方向取整函數(shù)，rand為（0，1）之間的隨機(jī)數(shù)。

Step10：根據(jù)Metropolis準(zhǔn)則判斷是否接受新解。令Δw＝w1-w0，當(dāng)Δw＜0，則接受w1作為新的當(dāng)前解w0， 否則以概率 exp（-Δw／T）接受w1作為新的當(dāng)前解w0，num＝num+1。

Step11：判斷當(dāng)前迭代次數(shù)num是否小于當(dāng)前溫度下的最大迭代次數(shù)maxnum。當(dāng)num＜max-num時，返回Step9；否則，執(zhí)行Step12。

Step12：用新解替換種群3中的原解。更新當(dāng)前溫度T，執(zhí)行降溫操作T＝low?T0，l＝l+1，并令num＝0。

Step13：判斷是否滿足算法終止條件。當(dāng)gen＝maxgen時，算法終止。否則，gen＝gen+1，返回Step3。

以目標(biāo)函數(shù)和罰函數(shù)之和作為Step3中的適應(yīng)度函數(shù)，如式 （28）所示：

3 仿真實(shí)驗(yàn)及分析

3.1 算例參數(shù)取值

現(xiàn)設(shè)計由1個主機(jī)廠、12個備選配送中心及40個4s店構(gòu)成的物流網(wǎng)絡(luò)仿真算例，由于篇幅限制僅列出6個備選配送中心及20個4s店的數(shù)據(jù)，設(shè)置備選配送中心相關(guān)參數(shù)如表1所示，備選配送中心與4s店的距離如表2所示，其他各項參數(shù)如表3所示，設(shè)置3種典型情景1、2、3分別代表需求量低、中、高的情況，需求均值及方差如表4所示。

表1 備選配送中心相關(guān)參數(shù)

表2 備選配送中心和4s店的距離 （km）

表3 其他參數(shù)取值

表4 不同情景下4s店需求均值及標(biāo)準(zhǔn)差數(shù)據(jù)

續(xù) 表

表5 魯棒優(yōu)化模型與隨機(jī)規(guī)劃模型性能對比

3.2 模型有效性分析

為了比較魯棒優(yōu)化模型B和隨機(jī)規(guī)劃模型C的性能，以上述數(shù)據(jù)為基礎(chǔ)設(shè)置3組不同節(jié)點(diǎn)規(guī)模的物流網(wǎng)絡(luò)，分別用模型B和模型C進(jìn)行求解。令情景g下的最優(yōu)目標(biāo)函數(shù)值為，魯棒優(yōu)化解在情景g下的目標(biāo)函數(shù)值為，魯棒優(yōu)化解的總成本為Zr，的相對遺憾值，ωr的方差為Vr。 隨機(jī)規(guī)化解在情景g下的目標(biāo)函數(shù)值為，隨機(jī)規(guī)劃解的總成本為Zs，的相對遺憾值，ωs的方差為Vs。Zr與Zs的相對差值為Zr-s＝（Zr-Zs）／Zs×100%，Vs與Vr的相對差值為Vs-r＝（Vs-Vr）／Vr×100%。 取p＝0.12， 結(jié)果如表5所示。

由表5可知，魯棒優(yōu)化解的相對遺憾值方差遠(yuǎn)小于隨機(jī)規(guī)劃解的相對遺憾值方差，兩者的相對差值最大達(dá)到了533.33%，而方差越小說明遺憾值的波動越小，在各種情景下的表現(xiàn)越穩(wěn)定，物流網(wǎng)絡(luò)的魯棒性越強(qiáng)［15］。且ωr的變動范圍為（0，11.48%），說明即使最壞的情景發(fā)生，魯棒優(yōu)化解在該種情景下的遺憾值也能控制在p值內(nèi)。而ωs的變動范圍為（0.03%，27.31%），在不同情景下的表現(xiàn)存在較大差距，這表明當(dāng)某種特定情景發(fā)生時 （如6?8算例中的情景1），隨機(jī)規(guī)劃解會產(chǎn)生較大的遺憾值，極大地增加了決策風(fēng)險。此外，雖然隨機(jī)規(guī)化模型的總成本較魯棒優(yōu)化模型低，但Zr-s最大為4.55%，成本差距并不明顯。因此，相比隨機(jī)規(guī)劃模型，魯棒優(yōu)化模型所求得的物流網(wǎng)絡(luò)優(yōu)化方案具有更強(qiáng)的魯棒性，能夠有效減少決策風(fēng)險。

3.3 算法有效性分析

為了驗(yàn)證本文所設(shè)計遺傳退火算法 （GASA）的有效性，以上述數(shù)據(jù)為基礎(chǔ)，設(shè)計不同規(guī)模算例，對比GASA與遺傳算法 （GA）及模擬退火算法 （SA）的求解結(jié)果，經(jīng)過調(diào)參實(shí)驗(yàn)后，設(shè)置遺傳算法的種群規(guī)模popsize為100，交叉概率Pc為0.9，變異概率Pm為0.06；模擬退火算法的初溫T0為9000，降溫系數(shù)low為0.99，每一溫度下的最大迭代次數(shù)maxnum為400。對于規(guī)模較小的算例1、2、3，遺傳算法和模擬退火算法的最大迭代次數(shù)均設(shè)置為1000，對于規(guī)模較大的算例4，設(shè)置最大迭代次數(shù)為5000。設(shè)置遺傳退火算法的相關(guān)參數(shù)取值與上述兩種算法保持一致，?＝0.05，θ＝1000000，在MATLAB R2017b環(huán)境下測試，每一算例均運(yùn)行10次，3種算法在情景1的運(yùn)行結(jié)果對比如表6所示。

表6 算法運(yùn)行結(jié)果對比

從表6可以看出，GASA在10次中的最優(yōu)值、最差值及平均值均小于SA與GA，雖然在融合了兩種算法的基本操作步驟后，GASA的平均運(yùn)行時間大于SA與GA，但能夠滿足物流網(wǎng)絡(luò)優(yōu)化問題對時效性的要求。通過以上的對比分析可證明本文所設(shè)計遺傳退火算法的有效性。

3.4 參數(shù)敏感性分析

保持其他參數(shù)取值不變，在6?10規(guī)模網(wǎng)絡(luò)（即i1～i6、j1～j10） 下對遺憾值限定系數(shù)p及碳稅φ分別進(jìn)行敏感性分析。

（1）遺憾值限定系數(shù)p敏感性分析

p值對目標(biāo)函數(shù)值的影響如圖3所示。由圖3可知，隨著p值的減少，目標(biāo)函數(shù)值逐漸增大，這意味著物流網(wǎng)絡(luò)的魯棒性與成本之間存在效益背反規(guī)律，若要增強(qiáng)網(wǎng)絡(luò)的魯棒性則需付出更高的成本。當(dāng)p＜0.1430時，算法未能找到可行解，這是因?yàn)閜存在一個下限值pmin，且當(dāng)p接近pmin時，目標(biāo)函數(shù)值對p的變化較為敏感，此時增強(qiáng)網(wǎng)絡(luò)的魯棒性所付出的成本相對更高。

而p＞0.1678時，目標(biāo)函數(shù)值不再發(fā)生變化，這是因?yàn)閜存在一個上限值pmax，當(dāng)p＞pmax時，約束條件 （24）將失效，此時即等同于隨機(jī)規(guī)劃模型。因此，決策者需要在魯棒性和成本之間均衡選擇，通過動態(tài)調(diào)整參數(shù)p調(diào)節(jié)整個物流網(wǎng)絡(luò)的魯棒性水平。

圖3 遺憾值限定系數(shù)敏感性分析

（2）碳稅值φ敏感性分析

調(diào)節(jié)碳稅φ的大小，得到碳排放量和總成本隨碳稅φ的變化情況如圖4所示。從圖4可以看出，隨著碳稅值的增加，碳排放量呈現(xiàn)階梯型遞減特征，目標(biāo)函數(shù)值則呈階段性線性增長趨勢——當(dāng)碳稅值在某一區(qū)間 （如φ＝0～0.11，0.12～0.20，0.21～0.30時）內(nèi)增加時，魯棒優(yōu)化方案及碳排放量均不發(fā)生變化，在該區(qū)間內(nèi)，目標(biāo)函數(shù)值呈線性增長，其斜率即對應(yīng)魯棒優(yōu)化方案的碳排放量。但當(dāng)碳稅值增加到某一臨界點(diǎn) （φ＝0.12，0.21）時，碳排放量突然下降。這是由于此時改變優(yōu)化方案所減少的碳稅成本大于所增加的其他項成本，企業(yè)將選擇碳排放量更少的優(yōu)化方案 （如表7所示），目標(biāo)函數(shù)值在該階段的斜率也隨之減小。因此，碳稅政策制定者需要根據(jù)這一規(guī)律，結(jié)合經(jīng)濟(jì)發(fā)展情況進(jìn)行相應(yīng)的調(diào)研分析，找到碳排放減少的臨界點(diǎn)設(shè)置碳稅值，使其在有效減少碳排放量的同時盡可能的降低企業(yè)的經(jīng)濟(jì)負(fù)擔(dān)。

圖4 碳稅值敏感性分析

表7 不同碳稅值下的優(yōu)化方案

4 結(jié) 論

根據(jù)以上分析，本文得出如下結(jié)論：

（1）建立了考慮碳稅和多種需求量情景的整車物流網(wǎng)絡(luò)魯棒優(yōu)化模型，通過與隨機(jī)規(guī)劃模型對比，證明魯棒優(yōu)化模型能夠有效減少決策風(fēng)險，但提高網(wǎng)絡(luò)的魯棒性必須付出更高的成本。

（2）設(shè)計了遺傳退火算法對模型進(jìn)行求解，并通過不同規(guī)模的算例實(shí)驗(yàn)證明其較參數(shù)取值相同下的遺傳算法和模擬退火算法性能更優(yōu)。

（3）靈敏度分析實(shí)驗(yàn)表明，總成本隨著碳稅值的增加逐漸增加，碳排放量隨著碳稅值的增加呈現(xiàn)明顯的階梯狀下降特征，碳稅政策制定者可利用這一規(guī)律在碳排放量變化的臨界點(diǎn)設(shè)置碳稅值，使其既能達(dá)到減排目的又能最大限度地減輕企業(yè)的經(jīng)濟(jì)負(fù)擔(dān)。