李魁梅，鄭 波

(1.重慶郵電大學(xué)移通學(xué)院 淬煉商學(xué)院，重慶 401520；2.重慶廣播電視大學(xué) 管理學(xué)院，重慶 401520)

隨著“一帶一路”國家戰(zhàn)略和“長江經(jīng)濟帶建設(shè)”的推進，經(jīng)濟和貿(mào)易全球化進程加快，對運輸?shù)囊笠灿鷣碛?，多式?lián)運在現(xiàn)代物流運輸中的地位越來越重要，成為當前研究的重點和熱點，取得了豐富的研究成果。分析相關(guān)文獻可知，相關(guān)研究主要集中在影響多式聯(lián)運路徑優(yōu)化的因素和路徑優(yōu)化算法研究兩方面。影響因素方面，隨著環(huán)保意識的增強，將尾氣排放、噪音污染等造成的運輸負外部性損失進行核算。Ricci等[1]將污染造成的外部成本考慮到運輸成本中，構(gòu)建了貨物運輸?shù)纳鐣杀居嬎隳Ｐ?，有重要參考價值。Macharis[2]提出運輸價格有效定價的途徑是外部成本內(nèi)部化，并在外部成本內(nèi)化這一前提下，研究了燃油價格波動對多式聯(lián)運服務(wù)范圍的影響。而國內(nèi)一般僅將二氧化碳排放作為外部性的考慮因素，如文獻[3]在進行多式聯(lián)運綜合滿意度建模時，將政府滿意度直接與二氧化碳排放量掛鉤。在模型方面，學(xué)者們分別提出了基于成本最小、運輸時間最少和運輸距離最短的多式聯(lián)運數(shù)學(xué)模型[4-6]。在算法方面，現(xiàn)有研究主要聚焦于智能算法，如文獻[7]用字符對個體進行編碼，提出了一種求解帶時間窗多式聯(lián)運問題的和聲算法。文獻[8]提出了遺傳算法和蟻群算法相結(jié)合的混合算法。文獻[9]則提出了基于固定優(yōu)先權(quán)編碼的遺傳算法，求得了長大貨物多式聯(lián)運問題的最優(yōu)解。

根據(jù)文獻研究發(fā)現(xiàn)，在進行多式聯(lián)運路徑優(yōu)化研究時，影響因素方面僅僅將時間因素作為約束條件，并沒有轉(zhuǎn)化入目標函數(shù)中，更沒有考慮尾氣排放外部成本(或者單純考慮二氧化碳排放)。此外，因為市場競爭激烈，越來越多的產(chǎn)品對時間敏感，比如筆記本電腦、時裝等，所以在進行路徑優(yōu)化時有必要增加考慮產(chǎn)品的價值特性，即衰變成本?；诖耍疚奶岢鰧⑦\輸過程中的排放、在途庫存、不準時(提前或逾期)和衰變構(gòu)建成本模型，結(jié)合運輸工具成本，獲得綜合運輸成本模型，進而以綜合運輸成本最低為目標構(gòu)建多式聯(lián)運路徑優(yōu)化模型。針對模型求解，因與一般智能算法相比，BA算法全局尋優(yōu)能力較強、穩(wěn)定性較好，不足之處在于收斂速度較慢，所以本文引入混沌機制、動態(tài)自適應(yīng)慣性權(quán)重、非均勻變異策略和反饋機制，獲得改進的HBA算法求解模型。

1 問題描述及數(shù)學(xué)模型

1.1 問題描述及假設(shè)

一批貨要從起點城市運輸?shù)浇K點城市，且兩者之間有若干城市可以中轉(zhuǎn)，且任意相連的兩個城市可能有多種運輸方式供選擇，除起點和終點外，每個節(jié)點均可實現(xiàn)貨物在不同運輸方式之間進行轉(zhuǎn)換。根據(jù)文獻[3]和[8]，結(jié)合各運輸方式市場供給情況，遵循合理假設(shè)原則，在如下假設(shè)條件(約束條件)下，合理規(guī)劃運輸方案，使綜合運輸成本最低。

1) 運輸和中轉(zhuǎn)工具運載能力充足，滿足貨物運輸和中轉(zhuǎn)需求；

2) 只考慮運輸?shù)某杀竞蜁r間，不考慮中轉(zhuǎn)的成本和時間；

3) 城市之間可以根據(jù)運輸條件自由選擇多種運輸方式，但每種運輸方式均不能超過該方式最大載重約束。

1.2 綜合運輸成本模型

綜合運輸成本由運輸會產(chǎn)生排放成本、在途成本、不準時成本、衰變成本和工具成本構(gòu)成。

1.2.1 排放成本

運輸中產(chǎn)生的排放物主要包括CO2、CO、NOX、PM、SO2、HC等6種，且不同運輸方式下各排放物的單位排放量和不同排放物的單位排放成本不一樣[10]。第k種運輸方式下，運輸單位貨物單位距離的排放成本如式(1)所示。

式中，k∈{1，2，···，K}，且K=3，分別表示公路、鐵路和水路運輸； ?∈{1，2，· ··，6}，分別表示CO2、CO、NOX、PM、SO2、HC 6 種排放物；Ek?表示第k種運輸方式第 ?種排放物單位排放量；Ce?表示第? 種排放物單位排放成本。

1.2.2 在途成本

貨物運輸過程中會占用大量資金，且產(chǎn)生保險、倉儲、稅費、搬運等費用，進而造成在途成本。根據(jù)文獻[11]，單位貨物單位時間產(chǎn)生的在途成本為

式中,GV表示貨物價值；H表示貨物單位時間持有成本。

1.2.3 不準時成本

不準時(提前或逾期)成本指未按照客戶規(guī)定時間送達貨物而產(chǎn)生的懲罰性費用為

式中，[ET, LT ]為客戶規(guī)定貨物送達的時間窗；tDC為貨物送達目的城市DC的時間；PE、PL分別為貨物早送達和晚送達的單位時間懲罰成本。

1.2.4 衰變成本

貨物價值通常會隨著時間衰減，進而產(chǎn)生衰變成本。文獻[11]用λ表示單位時間衰變率，則單位貨物單位時間衰變成本為

1.2.5 工具成本

工具成本指運輸過程中產(chǎn)生的運輸工具租賃費用、油耗費用和司機薪酬等。如果已知利用第k種運輸方式運輸單位貨物單位距離的工具成本為則再乘以運輸距離和貨物重量即可求得工具成本。

1.3 模型建立

I為運輸貨物總量；V為城市集合；m為城市節(jié)點數(shù)；dijk為城市i和j之間第k種運輸方式的運輸距離；tijk為城市i和j之間第k種運輸方式的運輸時間；分別為城市i和j之間第k種運輸方式的最大載重和實際載重；eijk為0-1變量，表示是否采用第k種運輸方式。

建立基于綜合運輸成本最小的多式聯(lián)運數(shù)學(xué)模型。

式(5)表示運輸成本最小的目標函數(shù)；式(6)表示起點城市運輸貨物量為輸入流量；式(7)表示所有貨物送達終點城市；式(8)表示每個城市的輸入流量與輸出流量相等；式(9)表示兩個城市之間任一運輸方式運輸量均不能超過該方式最大載重約束；式(10)表示城市間可選擇多種運輸方式；式(11)表示0-1變量約束。

2 基于混合蝙蝠算法的模型求解

2.1 標準蝙蝠算法

蝙蝠算法(BA)是一種基于迭代的優(yōu)化技術(shù)，初始化為一組隨機解，然后通過迭代搜尋最優(yōu)解，且在最優(yōu)解周圍通過隨機飛行產(chǎn)生局部新解，加強了局部搜索，是一種搜索全局最優(yōu)解的有效方法。與其他算法相比，蝙蝠算法(BA)在準確性和有效性方面遠優(yōu)于其他算法。蝙蝠算法(BA)包括搜索和尋優(yōu)兩大機制[12-13]。

2.1.1 搜索機制

搜索機制主要利用頻率h對蝙蝠位置和速度進行更新。在D維空間中，用和(xb1,xb2,···,xbD)分別表示蝙蝠b在t時刻的速度和位置，則t+1時刻更新為

式中，X?表示當前最優(yōu)蝙蝠所在位置；hb=(hmax?hmin)φ(φ為服從均勻分布的隨機數(shù)，取值范圍為[0,1])。

2.1.2 尋優(yōu)機制

通過尋優(yōu)機制，每只蝙蝠位置Xold在最優(yōu)個體附近搜索，并根據(jù)式(14)更新位置。

2.2 編碼及解碼方式

2.2.1 全局流量按比例分配

式中，σi表示第i條路徑的貨物流量比例；xi表示第i條路徑的貨物流量；C1s表示與第i條路徑同一等級的所有路徑集合(輸出節(jié)點為s)。則第i條路徑承擔(dān)的貨物流量為

式中，Oflows表示路徑i對應(yīng)輸出節(jié)點s的輸出流量值。如圖1所示，路徑4對應(yīng)流量輸出節(jié)點1，且同一等級路徑集為C11={1,3,4,5,6}。 因此，σ4=x4/用Oflow1表示節(jié)點1輸出流量，則路徑4的流量全局流量的按比例分配的具體流程如圖2所示。

圖1 全局流量按比例分配Figure 1 Proportional distribution of global flow

圖2 全局流量分配流程Figure 2 Distribution process of global flow

2.2.2 局部流量返上級分配

實際上，每條路徑的貨物承載能力有限，簡單按比例進行流量分配，必然產(chǎn)生較多不可行方案，進而影響算法求解效果。對此，本文設(shè)計了將節(jié)點超額流量返回上級節(jié)點進行再次分配的思路。具體操作流程為：判斷節(jié)點輸入流量和最大流量約束大小，如果輸入流量大于最大流量約束，則需要進行局部調(diào)整計算節(jié)點超額流量，即輸入流量和最大流量約束之差；將超額流量返回上一節(jié)點，按滿載率對其他可行路徑進行排序，依次分配超額流量，優(yōu)先滿足單條路徑滿載。根據(jù)前文可知，節(jié)點4 輸入流量為f4，其對應(yīng)輸出路徑為7。用flimit,i表示路徑i的最大流量約束，那么如果f4＞flimit,7，則節(jié)點4流量超額，需要進行局部調(diào)整。如圖3所示，用fleft4表示節(jié)點4的超額流量，取值為f4?flimit,7。將fleft4返回節(jié)點1進行再次分配，可知節(jié)點1對應(yīng)的輸出路徑還有1、3、5和6。計算每條路徑滿載度 De=flimit?f，并進行降序排序，然后將fleft4進行依次分配。用qi(i∈{1,3,5,6})表 示各條路徑再次分配流量，則0≤qi≤Dei,如圖4所示。通過上述過程，便能實現(xiàn)超額流量的再次分配，大量減少不可行方案的產(chǎn)生。

圖3 超額流量返回上級節(jié)點Figure 3 Excess flow returned to parent node

圖4 上級節(jié)點超額流量再分配Figure 4 Reallocation of excess flow in superior nodes

2.3 設(shè)計適應(yīng)度函數(shù)

適應(yīng)度值是評判個體優(yōu)劣的重要標準，本文根據(jù)式(5)設(shè)計算法適應(yīng)度函數(shù)為

式中，η為正值常數(shù)。

2.4 引入混沌狀態(tài)

智能算法初始種群質(zhì)量決定著算法的收斂速度和求解精度[6]。標準BA算法在缺乏有關(guān)先驗信息的情況下，采用隨機方式生成初始種群，無法確保初始種群在解空間中的均勻分布。對此，本文利用混沌系統(tǒng)的隨機性、遍歷性和規(guī)律性來優(yōu)化蝙蝠初始位置，提高初始種群質(zhì)量，增強算法求解性能。具體包括兩個步驟：首先，根據(jù)式(20)構(gòu)造大量蝙蝠位置和速度向量；然后擇優(yōu)選擇部分個體組成蝙蝠初始種群。

式(20)為Logistics映射(迭代)，μ為控制參數(shù)。當 μ=4時 ，y0∈(0,1)，表示Logistics處于完全混沌狀態(tài)。

2.5 引入動態(tài)自適應(yīng)慣性權(quán)重

根據(jù)式(12)可知，標準BA 算法速度更新系數(shù)恒定為1，不但大大降低了種群多樣性和個體靈活性，還導(dǎo)致算法在尋優(yōu)過程中不能很好地平衡局部搜索和全局搜索的關(guān)系，進而影響算法尋優(yōu)效率。針對該問題，本文提出一種動態(tài)自適應(yīng)慣性權(quán)重系數(shù)，具體表達式如式(21)所示。

式中，ωmax、 ωmin分別為 ω的最大值和最小值；φ表示服從均勻分布的隨機數(shù)，取值范圍為[0,1]；τ表示慣性偏離因子，取值范圍為[0.1,0.9]； β表示服從貝塔分布的隨機數(shù)，取值為范圍為(0,1)。分析式(21)( ωmax?ωmin)φ使 ω能 夠 在[ωmin, ωmax]內(nèi) 隨 機 取值，而貝塔分布隨機數(shù) β讓ω 隨 機偏離，使得 ω的取值分布更加均勻和靈活。同時，在 β前加入τ，很好地控制了 ω的偏離幅度?？梢钥闯?，引入動態(tài)自適應(yīng)權(quán)重，不但可以增加代與代之間的蝙蝠隨機性，也增加了每一代內(nèi)蝙蝠的隨機性，提高了蝙蝠種群多樣性。

在標準BA算法中引入動態(tài)自適應(yīng)慣性權(quán)重，則算法速度更新公式變?yōu)?/p>

2.6 引入非均勻變異策略

針對標準BA算法容易陷入局部最優(yōu)解的問題，本文采用非均勻變異操作對其進行改進。與遺傳算法變異算子不同，非均勻變異操作能夠根據(jù)迭代次數(shù)自適應(yīng)調(diào)整搜索步長，進而使算法在進化過程中始終保持跳出局部最優(yōu)解的能力[13]。如式(23)所示，如果對蝙蝠個體Xi=(xi1,xi2,···,xin,···,xiD)進行變異操作，則第n維位置向量變異為

式中，UB、LB分別表示當前迭代中所有蝙蝠個體第n維變量的最大值和最小值；(z,y)=y(1?φ(1?z/Z)Λ)，表示非均勻變異步長，z表示當前迭代次數(shù)，Z表示最大迭代次數(shù)；φ表示服從均勻分布的隨機數(shù)，取值范圍為[0,1]；Λ表示非均勻度參數(shù)。分析可知，通過非均勻變異操作，算法在迭代初期擁有較大的搜索范圍，提高了全局尋優(yōu)能力。隨著迭代次數(shù)的增加，搜索半徑逐漸減小，算法只能在當前解附近進行搜索，提高了算法搜索精度。

2.7 引入反饋機制

在標準BA 算法中，蝙蝠個體通過不斷更新自身位置和速度來靠近獵物，運動方式比較單一[10]。對此，本文引入信息反饋機制，即每次迭代中最優(yōu)個體與最差個體進行信息反饋，通過這種信息交流，引導(dǎo)距離獵物較遠的個體快速抵達獵物附近，進而提高算法收斂速度。信息反饋機制如式(24)所示。

式中，Xworst、X?分別表示當前最差和最優(yōu)蝙蝠所在位置；?表示服從均勻分布的隨機數(shù)，取值范圍為[0,1]。

2.8 模型求解步驟

根據(jù)前文所述，采用混合蝙蝠算法HBA求解多式聯(lián)運問題的流程如下。

步驟1設(shè)置算法參數(shù)，包括候選種群規(guī)模M、初始種群規(guī)模N、最大迭代次數(shù)Z、最大和最小慣性權(quán)重ωmax、ωmin等。

步驟2隨機生成一個實數(shù)數(shù)組 (x1,s1)并按照式(20)生成其他D-1維向量，作為蝙蝠b的位置和速度向量 (Xb,Sb)。同理，可生成其他M-1個蝙蝠個體。對于Xb，根據(jù)3.2節(jié)進行解碼，并根據(jù)式(19)測算適應(yīng)度值C(Xb)。擇優(yōu)選擇N個蝙蝠個體組成初始種群。

步驟3記錄當前種群最優(yōu)個體X?及目標函數(shù)值gC。

步驟4判斷是否迭代了Z次。是，輸出種群最優(yōu)個體X?，算法停止；反之，進入步驟5。

步驟5根據(jù)式(21)和(22)更新蝙蝠個體的速度向量，再按照式(13)更新位置向量。

步驟6生成隨機數(shù)rand1，滿足條件，則按照式(23)更新蝙蝠個體i得到Xnew，進入步驟7；反之，進入步驟8。

步驟7生成隨機數(shù)rand2，同時滿足rand1>A(b)和C(Xnew)

步驟8按照式(15)和(16)更新A(b)和r(b)。

步驟9選擇種群最差個體，按照式(24)進行更新，返回步驟4。

3 仿真分析

本文采用包括公路、鐵路和水路3種運輸方式的隨機運輸網(wǎng)絡(luò)來驗證HBA算法的有效性。如圖5所示，該運輸網(wǎng)絡(luò)共有10個節(jié)點，節(jié)點1為起點城市，節(jié)點10為終點城市。不同節(jié)點之間的運輸方式可以分為單一運輸和多式聯(lián)運，如節(jié)點4到節(jié)點6只有鐵路運輸，節(jié)點7到節(jié)點8包括公路運輸和鐵路運輸，可以同時選擇2種方式進行運輸。運輸網(wǎng)絡(luò)中，括號內(nèi)數(shù)值分別表示不同運輸方式的運輸距離(單位為km)和最大容量限制(單位為t)。一般來說，公路運輸運載能力不存在約束。

圖5 公路、鐵路、水路聯(lián)運網(wǎng)絡(luò)Figure 5 Intermodal transport network of Road, railway and waterway

企業(yè)需要從節(jié)點1運輸重量500 t、價值880萬元的貨物到節(jié)點10，節(jié)點10的時間窗為[18,23]，不準時懲罰成本為每小時0.45 萬元，單位貨物持有成本系數(shù)為每小時0.0003%，單位貨物衰減率為每小時0.018%，公路、鐵路、水路運輸速度分別為80 km/h、100 km/h和54 km/h，運輸工具成本分別為80元/( t·km)、40元/( t·km)和16元/ (t·km)，不同運輸方式排放量和排放成本[14]如表1所示。

表1 3種運輸方式平均排放量和排放成本Table 1 Average emissions and emission costs of three transportation modes

利用HBA算法對上述算例進行50次求解，收斂情況如圖6所示，最優(yōu)運輸方案如圖7所示(圖中數(shù)值表示對應(yīng)運輸方式運輸貨物的重量),最優(yōu)運輸成本為2.56萬元。由此可見，多式聯(lián)運路徑優(yōu)化模型可行有效且HBA算法能夠適用于多式聯(lián)運問題的求解。

為進一步驗證HBA算法的求解性能，分別采用HSA算法[7]和HA算法[8]對上述算例進行50次求解，最優(yōu)值對應(yīng)收斂狀況如圖8所示，仿真結(jié)果如表2所示。

圖6 HBA算法收斂圖Figure 6 Convergence graph of HBA algorithm

圖7 最優(yōu)運輸方案Figure 7 The optimal transportation scheme

圖8 不同算法收斂對比Figure 8 Convergence comparison of different algorithms

表2 不同算法仿真結(jié)果對比Table 2 Comparison of simulation results between different algorithms

分析圖8可知，與HSA算法和HA算法相比，HBA算法收斂較快，且求得運輸成本最低，說明HBA算法具有較好的求解性能。

由表2可知，在全局尋優(yōu)能力上，HBA算法能夠求得最優(yōu)值，且最優(yōu)值搜索成功率高達90%，HSA算法和HA算法最優(yōu)值搜索成功率均為0，故HBA算法全局尋優(yōu)能力最強。在算法穩(wěn)定性上，HBA算法求得最劣值、平均值依次優(yōu)于HA算法和HSA算法，所以算法穩(wěn)定性由強到弱排序為HBA算法、HA算法和HSA算法。在算法運行速度上，各算法平均迭代次數(shù)和計算時間從小到大依次為HBA算法、HSA算法和HA算法，可見HBA算法在提高算法全局尋優(yōu)能力和穩(wěn)定性的同時，保持了較快的運行速度。綜上所述，HBA算法在求解多式聯(lián)運問題上優(yōu)于HSA 算法和HA算法。

4 結(jié)論

1) 為更加客觀全面反映運輸生產(chǎn)實際，將客戶滿意度和企業(yè)責(zé)任相結(jié)合，構(gòu)建考慮運輸工具尾氣排放、在途庫存損失、不準時損失和衰變損失以及傳統(tǒng)意義的運輸費用5個方面的多式聯(lián)運路徑優(yōu)化模型。

2) 模型計算采用混合蝙蝠算法(HBA)求解，算法在標準BA算法基礎(chǔ)上，改進了算法解碼方式，引入了混沌機制、動態(tài)自適應(yīng)慣性權(quán)重、非均勻變異策略和反饋機制。

3) 通過算例分析表明，多式聯(lián)運路徑優(yōu)化模型可行有效，且HBA算法在全局尋優(yōu)能力、算法穩(wěn)定性和運行速度上均優(yōu)于HSA算法和HA算法，有效提高了蝙蝠算法求解多式聯(lián)運問題的性能。