李子慷, 劉林冬, 于成成

(中國科學(xué)技術(shù)大學(xué) 管理學(xué)院國際金融研究院,安徽 合肥 230026)

0 引言

最基本的設(shè)施選址問題分為無容量限制和有容量限制兩種。其中無容量限制設(shè)施選址問題(Uncapacitated Facility Location Problem, UFLP)在文獻中得到了廣泛的討論[1]。Krarup對研究UFLP問題的文獻進行了綜述，并指出了UFLP問題是NP-難的[2]。逆優(yōu)化問題最初由Tarantola給出了逆優(yōu)化問題在地球物理科學(xué)中的詳細討論[3]。Heuberger給出了關(guān)于逆優(yōu)化問題較為全面的綜述[4]。Ahuja和Orlin證明了在范數(shù)L1和范數(shù)L∞下，如果原線性規(guī)劃問題多項式可解則其逆問題也是多項式可解的[5]。通常來講，只有逆問題的可行域是一個多面體，才能確保在范數(shù)L1下的逆問題在多項式時間內(nèi)是可解的[6]。Schaefer給出了一般整數(shù)規(guī)劃逆問題的多面體描述，并提出了兩種算法求解[7]。

由于設(shè)施選址逆問題在一般網(wǎng)絡(luò)下，大多數(shù)是NP-難的，很多學(xué)者考慮了特殊情況下的逆選址問題。Cai等人考慮了逆單中心選址問題，證明了即使在單中心選址原問題是多項式可解的情況下，逆單中心選址問題仍然是NP-難的[6]。Alizadeh等人考慮了在樹狀圖中邊長度僅允許增加時的逆單中心選址問題，給出了一個時間復(fù)雜度為O(nlogn)的精確算法[8]。Guan給出了在加權(quán)L∞范數(shù)下，逆1-中值選址問題線性時間算法，并證明了在加權(quán)哈明距離下該逆問題是NP-難的[9]。

在逆選址值問題上，Berman等人證明了在二分網(wǎng)絡(luò)中的逆選址值問題是NP完全的，同時給出了在樹形網(wǎng)絡(luò)中該問題的線性表達式[10]。Zhang等人研究樹形網(wǎng)絡(luò)中的逆選址值問題，并給出了一個時間復(fù)雜度為O(n2logn)的多項式時間算法，同時提出該問題的變體，證明了逆多設(shè)施選址值問題可以轉(zhuǎn)化為幾個逆單設(shè)施選址值問題[11]。

設(shè)施選址逆問題具有潛在的應(yīng)用前景，它為我們展示了如何花費盡可能少的成本，從而固定開啟某些現(xiàn)有設(shè)施以提供便利。對應(yīng)到模型中，即給定了開設(shè)哪些設(shè)施，這些設(shè)施對哪些顧客進行服務(wù)，目標(biāo)函數(shù)則為使得設(shè)施和運輸成本改變量達到最小。因而研究最基本的無容量限制設(shè)施選址逆問題及其求解算法是十分具有現(xiàn)實價值的。

本文研究了無容量限制設(shè)施選址逆問題(Inverse Uncapacitated Facility Location Problem，IUFLP)模型的建立與求解方法。用行生成算法將IUFLP問題分為主問題和子問題，其中子問題為一般的無容量限制設(shè)施選址問題(UFLP)。對子問題進行不同的處理方式，可以分別得到IUFLP問題的上下界。

本文安排如下：第一節(jié)給出了UFLP的定義以及IUFLP模型的建立。第二節(jié)給出了IUFLP的行生成求解方法。第三節(jié)對由行生成算法得到的子問題進行了改進，得到能求出該逆問題上下界的啟發(fā)式算法。第四節(jié)給出了IUFLP的簡單示例以及計算結(jié)果。第五節(jié)給出了對于該逆問題延拓以及算法改進的展望。

1 UFLP和IUFLP問題的模型

設(shè)施選址問題通常有如下的定義：假設(shè)有m個設(shè)施和n個顧客。我們希望選擇(1)哪些設(shè)施要打開，以及(2)哪些打開的設(shè)施要用于向哪些顧客提供服務(wù)，以便以最低的成本滿足需求。引入以下記號，令fi表示開啟設(shè)施的固定成本，i∈M,M={1,…,m}。rij表示運輸商品從設(shè)施i到顧客j的成本，j∈N,N={1,…,n}。無容量限制指每個設(shè)施可以提供的容量K是充分大的。

設(shè)為fi維m的行向量，將rij按行首尾相接為一個(m×n)維的行向量，則可用一個(m+m×n)維行向量c=(f,r)表示成本。

則無容量限制設(shè)施選址問題(UFLP)定義如下:

uij,vi∈{0,1},?i∈M,?j∈M

用更為簡潔的方式改寫包含足夠大容量K的上式，得到下式:

uij≤vi,?i∈M,?j∈N

uij,vi∈{0,1},?i∈M,?j∈N

在上述整數(shù)規(guī)劃模型中，uij表示顧客j是否被設(shè)施i所服務(wù)，如果被服務(wù)，則uij=1，否則uij=0；vi表示設(shè)施i是否開啟，如果開啟，則vi=1，否則vi=0。(1)式為UFLP的目標(biāo)函數(shù)，要使得開啟設(shè)施和運輸商品的總成本最??；(2)式表示每個顧客j∈N必須被一個設(shè)施所服務(wù)；(3)式表示如果顧客j被設(shè)施i所服務(wù)，則設(shè)施i必須開啟。

逆問題通常有如下的定義，當(dāng)給定一個加權(quán)優(yōu)化問題和一個固定可行解的實例，其相應(yīng)的逆問題就是盡可能少地(在一個固定范數(shù)下)修改加權(quán)參數(shù)，使得給定的解在該實例中成為原優(yōu)化問題的最優(yōu)解。

根據(jù)逆問題的定義，一般優(yōu)化問題的逆問題有如下形式：

mind‖d-c‖

s.t.d∈G(x0)

G(x0)={d∈Rn

|min{g(d,x)|∈D}=g(d,x0)

其中g(shù)(c,x)為原優(yōu)化問題的目標(biāo)函數(shù)，c∈Rn為參數(shù)向量，D為x的可行域，d∈Rn為調(diào)整后的參數(shù)向量。G(x0)表示使給定的可行解x0成為g(d,x0)最優(yōu)解的參數(shù)向量d的集合。

則L1范數(shù)下無容量限制設(shè)施選址的逆問題(IUFLP)可以定義如下:

minw×|d-c|T

s.t.dx0=U(d)

d∈Rm+mn

其中w=(wf,wr)為(m+mn)維行向量，d=(f0,r0)為調(diào)整成本后的(m+mn)維行向量，f0,r0分別表示調(diào)整之后的設(shè)施成本和運輸成本，c=(f,r)為上述的原成本參數(shù)，|d-c|T為成本調(diào)整量(絕對值)的(m+mn)維列向量，x0為給定原問題的一個可行解，U(d)為(1)式中參數(shù)為d時對應(yīng)的最小值。

由約束式可以看出，x0是滿足調(diào)整成本后UFLP所有約束的最優(yōu)解。對于每一個顧客j∈N可以被任一設(shè)施服務(wù)，共有|M|=m種選擇。因而對于所有顧客而言一共會有mn種不同情況。即滿足UFLP約束條件的可行解有mn個。同時定義滿足UFLP約束條件的mn個可行解構(gòu)成的解集為Y。因此從逆問題的定義出發(fā)建立求解模型，Y中包含可行解的個數(shù)是指數(shù)個。事實上，無容量限制設(shè)施選址的逆問題(IUFLP)是NP-難的。

定理1無容量限制設(shè)施選址的逆問題(IUFLP)是NP-難的。

2 行生成算法求解

從上述無容量限制設(shè)施選址逆問題的定義過程中，我們可以看出原優(yōu)化問題和其逆問題有著一定的聯(lián)系，即逆問題中的給定解對應(yīng)的成本需要小于等于所有滿足原問題約束條件的解對應(yīng)的成本。行生成算法作為一種有效求解規(guī)模較大線性規(guī)劃問題的方法。我們可以將行生成算法應(yīng)用于該問題中，并且容易發(fā)現(xiàn)此時的子問題即是UFLP原問題，因而行生成算法可以用來求解無容量限制設(shè)施選址的逆問題。

無容量限制設(shè)施選址的逆問題模型為:

?j∈N,(v′,u′)∈Y

其中Y同上述，是滿足原問題的所有可行解x′=(v′,u′)的集合。為方便顯示我們的結(jié)果，我們主要考慮逆最優(yōu)解問題，即逆問題。而對于逆最優(yōu)值問題，即給定V0而不是x0，根據(jù)式V0=cx0也可用相同思路求解。下文只針對逆問題進行分析求解。

將行生成算法應(yīng)用于IUFLP，將c對應(yīng)的設(shè)施和運輸部分用(f,r)分別表示。根據(jù)Ahuja在L1范數(shù)下逆問題的處理方法[5]，記c=(f,r)為原成本向量，d=(f0,r0)為調(diào)整后的成本向量。令

則ai,hij分別表示相對原成本的增加量，fi,rij表示相對原成本fi,rij的減少量。

可以得到IUFLP的主問題為:

?i∈M,?j∈N

ai≥0,bi≥0,?i∈M

hij≥0,lij≥0,?i∈M,?j∈N

(v′,u′)∈Ⅱ

其中(v′,u′)屬于一個可行解限制集Ⅱ?Y，限制集中的元素由子問題不斷生成加入其中。而主問題得到的(a′,b′,h′,l′)在每次迭代循環(huán)中帶入到子問題中求解。

子問題如下:

子問題求出得到(v′*,u′*)，如果目標(biāo)函數(shù)值小于0，將得到的(v′*,u′*)加入到主問題的限制集中。如果目標(biāo)函數(shù)值不小于0，程序結(jié)束。注意到子問題實際上就是UFLP。

求解IUFLP問題的行生成算法詳細步驟如下:

輸入：1)可行解：x0=(v0,u0)； 2)原設(shè)施成本：c=(f,r)；

輸出：1)最優(yōu)解：(a*,b*,h*,l*)；2)最小的總成本改變量：

3)對應(yīng)原成本的改變量：

初始化：初始限制集包含對于原問題的可行解即可。為方便且不再重新生成可行解，可取初始限制集為只含給定的可行解這一個元素的集合即I0={(v0,u0)}。

步驟1求解此時的主問題得到最優(yōu)解(a′,b′,h′l′)。

步驟2將步驟1得到的最優(yōu)解(a′,b′,h′l′)帶入到子問題中，求解子問題這一整數(shù)規(guī)劃得到最優(yōu)解x′*=(v′*,u′*)和目標(biāo)函數(shù)值T。

步驟3如果目標(biāo)函數(shù)值T<0，將得到的x′*加入到限制集中，Ⅱ=Ⅱ∪{x′*}回到步驟 1，重新求解主問題；如果目標(biāo)函數(shù)值T≥0，則結(jié)束算法，得到最終的限制集，記為ⅡI。根據(jù)此限制集ⅡI求解主問題，即可得到最優(yōu)解(a*,b*,h*,l*)。

3 啟發(fā)式算法得到上下界

由上述的行生成方法，我們可以知道要求得逆問題的最優(yōu)解就需要得到盡可能精確的主問題的限制集，即主問題中的約束條件。由于主問題的限制集是由子問題不斷生成得到的，可以看出行生成方法的核心是子問題的求解效率，即子問題求解越精確，逆問題的最優(yōu)解越精確。但由于子問題是設(shè)施選址的原問題，因而可以通過對子問題進行松弛，犧牲一定求解精度的情況下簡化求解過程。

當(dāng)對子問題進行如下幾種處理方法時，我們可以相應(yīng)得到行生成方法的幾點結(jié)論：

處理方法1：嚴(yán)格求解子問題(S)。

定理2當(dāng)嚴(yán)格求解子問題時，循環(huán)結(jié)束時此時得到IUFLP的最優(yōu)解。

說明：注意到從子問題中篩選出的解的數(shù)量依賴于給定的初始解，雖然行生成方法可以極大地減少原有的約束，但究竟需要篩選出多少解是無法判斷的。事實上，如果可以判斷出能夠篩選出多少解，根據(jù)篩選出的解構(gòu)成的解集直接求解主問題就可以解決IUFLP。

處理方法2：對子問題利用鄰域搜索方法得到啟發(fā)式整數(shù)解，同時設(shè)置一定的循環(huán)次數(shù)提前終止子問題的求解。

定理3對子問題利用鄰域搜索方法求解，并且設(shè)定一定的循環(huán)次數(shù)提前終止子問題的求解，此時求解得到IUFLP的下界。

對于利用鄰域搜索求解子問題的說明：鄰域搜索采用基本的本地搜索方法[12]，即對于當(dāng)前給定解采用增加，交換和減少的方式得到局部最優(yōu)整數(shù)解。

下面對鄰域搜索的三種方式進行說明。設(shè)此時未開啟設(shè)施的數(shù)量為p，開啟設(shè)施的數(shù)量為(m-p)。增加，即對于當(dāng)前未開啟的設(shè)施中開啟一個設(shè)施，作為一個增加鄰域，共有p個鄰域。交換，即對于當(dāng)前未開啟的設(shè)施中開啟一個設(shè)施，在開啟的設(shè)施中關(guān)閉一個設(shè)施，作為一個交換鄰域，共有p(m-p)個鄰域。減少，即對于當(dāng)前開啟的設(shè)施中關(guān)閉一個設(shè)施，作為一個減少鄰域，共有(m-p)個鄰域。

在一次求解過程中，對于三種方式得到的鄰域解對應(yīng)的成本進行比較，取最小值求得對應(yīng)的局部最優(yōu)解。將當(dāng)前給定解轉(zhuǎn)移到得到的局部最優(yōu)解，同時將該局部最優(yōu)解加入到主問題的限制集中，再進行下一次求解。

求解UFLP下界的啟發(fā)式算法詳細步驟如下：

輸入：1)可行解：x0=(v0,u0)；2)原設(shè)施成本：c=(f,r)；

輸出：1)最優(yōu)解：(a*,b*,h*,l*)；2)最小的總成本改變量：

3)相對原成本的改變量：

初始化：初始限制集包含對于原問題的可行解即可。為方便且不再重新生成可行解，可取初始限制集為給定的可行解即Ⅱ={(v0,u0)}。

步驟1求解此時的主問題得到最優(yōu)解(a′,b′,h′,l′)。

步驟2將步驟1得到的最優(yōu)解(a′,b′,h′,l′)帶入到子問題中，鄰域搜索的子問題目標(biāo)函數(shù)與前述保持一致，對給定的初始解采用增加，交換，減少操作得到局部最優(yōu)解x′*=(v′*,u′*)和目標(biāo)函數(shù)值T。

步驟3如果目標(biāo)函數(shù)值T<0，將得到的x′*加入到限制集中，Ⅱ=Ⅱ∪{x′*}回到步驟1，重新求解主問題；如果目標(biāo)函數(shù)值T≥0，則結(jié)束算法。此時得到最優(yōu)解(a*,b*,h*,l*)。

其中鄰域搜索中的增加，交換，減少操作如圖1所示：

圖1 鄰域搜索對當(dāng)前解進行三種操作的示意圖

處理方法3：對子問題(S)進行線性松弛求解。

定理4將子問題(S)進行線性松弛得到線性規(guī)劃問題，使用行生成算法求解則可以得到對應(yīng)逆問題的上界。

此時對應(yīng)的子問題模型變?yōu)?

行生成算法及使用不同處理方法處理子問題求解的流程示意圖如下:

圖2 行生成算法及不同處理方法流程示意圖

4 示例及計算結(jié)果

4.1 簡單的例子

我們舉一個簡單的例子來具體展示無容量設(shè)施選址逆問題。

圖3 原UFLP示例圖

圖4 給定可行解的IUFLP示例圖

如上圖所示，Ⅰ、Ⅱ、Ⅲ、Ⅳ為設(shè)施編號，每個設(shè)施開啟成本為10。a,c,b,d為顧客編號，設(shè)施與顧客有連線表示該顧客可以被相應(yīng)設(shè)施服務(wù)，連線上的數(shù)字表示對應(yīng)的運輸成本。黑色的細線表示該問題所有可行的服務(wù)連線，帶箭頭的粗線表示服務(wù)最優(yōu)解連線。相應(yīng)的設(shè)施開啟成本和運輸成本為:

圖3中給出了在給定設(shè)施開啟成本和運輸成本后的最優(yōu)解服務(wù)連線。圖4則是給定了如圖所示的服務(wù)連線，需要調(diào)整設(shè)施開啟成本或運輸成本使得該服務(wù)連線成為調(diào)整后最優(yōu)的服務(wù)連線。

能夠看出，原問題最優(yōu)解如圖3所示，即為開設(shè)設(shè)施Ⅱ服務(wù)顧客b，以及開設(shè)設(shè)施Ⅲ服務(wù)顧客a,c,d，對應(yīng)的最優(yōu)成本為27。而當(dāng)給定可行解如圖4所示時，即為開設(shè)設(shè)施Ⅲ服務(wù)顧客a,d，以及開設(shè)設(shè)施Ⅳ服務(wù)顧客b,c時，對應(yīng)的成本為29。需要改變方框中的運輸成本使得給定的可行解成為改變成本后的最優(yōu)解，即將方框中1改為2，方框中2改為3。此時最小改變量為2，可以使得給定的可行解成為改變成本后問題的最優(yōu)解。

4.2 計算結(jié)果

本文使用Windows 10操作系統(tǒng)，內(nèi)存8 GB處理器為Intel Core i7- 8700的PC電腦上進行了所有的數(shù)值實驗。所有算法均通過Matlab R2019a調(diào)用Gurobi求解器去實現(xiàn)。在數(shù)值實驗中，我們對于6種不同規(guī)模隨機生成了對應(yīng)規(guī)模下20種可能的實例。(m,n)=(10,10),(10,20),(10,30),(20,20),(20,40),(30,30)。

其中設(shè)施啟動成本和運輸成本均為隨機產(chǎn)生的整數(shù)值，且取值范圍分別為：

fi=[100,200]；rij=[1,100],[100,200]

在實際中往往不能改變設(shè)施啟動成本時，因而在接下來的計算中設(shè)置對于設(shè)施成本和運輸成本調(diào)整量的權(quán)重分別為wf=∞,wr=1。

我們將采用下列方式來比較最優(yōu)改變量和啟發(fā)式算法計算得到的改變量之間的差距。當(dāng)最優(yōu)改變量不為0時，我們采用相對值(上下界改變量-最優(yōu)改變量)/(最優(yōu)改變量) 的方式來比較gap的大小。相應(yīng)的結(jié)果見表1。由于給定的可行解是隨機生成的，因而給定可行解剛好是最優(yōu)解的可能性幾乎為零。而當(dāng)給定可行解剛好是最優(yōu)解時，根據(jù)上述的算法可以知道不管是精確求解算法還是啟發(fā)式算法，此時的最優(yōu)改變量均為0。但在實際問題中，幾乎不存在給定的可行解剛好是最優(yōu)解，因為這樣意味著不需要對原有成本做任何的改變，這里不再考慮這一極端情況。

結(jié)果如下:

表1 給定可行解非最優(yōu)解時,IUFLP問題的計算結(jié)果

在表中，(m,n)-1表示IUFLP問題的規(guī)模。例如(20,40)-1表示20個設(shè)施和40個顧客：-1表示fi=[100,200],rij=[1,100];-2表示fi=[100,200],rij=[100,200]。

OS表示嚴(yán)格求解子問題時得到的逆問題的最優(yōu)值。UB表示使用線性松弛的方法得到的 P的上界。LB表示使用鄰域搜索得到的 P的下界。在表1中GAP計算方法為相對比值法。GAP1表示UB與OS之間的gap,GAP2表示LB與OS之間的gap。T1，T2，T3分別為得到最優(yōu)值，上界，下界的計算耗時，單位為秒。表中所有的結(jié)果均是一個規(guī)模中產(chǎn)生的20個實例的平均值，并經(jīng)過四舍五入處理。

特殊說明：在使用處理方法二得到下界的過程中，設(shè)置了最大循環(huán)次數(shù)1000。而在處理方法一和三中，設(shè)置的最大循環(huán)次數(shù)較大不會提前終止程序。

計算結(jié)果表明，問題的規(guī)模(m,n)以及設(shè)施開設(shè)成本fi與運輸成本rij之間的大小關(guān)系不會對gap的大小造成太大影響。同時，由于該問題是NP-難的，隨著規(guī)模的變大，計算復(fù)雜度和計算時間都是指數(shù)增長的。根據(jù)計算結(jié)果可以看出，線性松弛得到的上界與最優(yōu)的改變量差距較小，但求解效率卻沒有較大提升。原因在于求解最優(yōu)的改變量時是使用整數(shù)規(guī)劃模型進行求解，而采用的求解器對于整數(shù)規(guī)劃的求解效率較高；同時線性松弛方法在每次迭代循環(huán)過程中產(chǎn)生的是非整數(shù)解，迭代循環(huán)次數(shù)相對變多也增大了求解需要的時間。而利用基本的本地搜索方法求出的下界與最優(yōu)的改變量非常接近，同時在求解效率上也有較大提升。因此，當(dāng)問題規(guī)模不大時，可以采用最優(yōu)解的方式進行求解；而當(dāng)問題規(guī)模較大時，可以采用啟發(fā)式方法結(jié)合一定的迭代循環(huán)次數(shù)得到不錯的下界。迭代循環(huán)次數(shù)越多，計算耗時就越多，但得到的下界就越好。在實踐中，針對具體的問題具體分析，使用啟發(fā)式方法同時設(shè)置一定的迭代循環(huán)次數(shù)可以得到不錯的成本改變量。

5 結(jié)論

本文針對無容量限制的設(shè)施選址逆問題進行了建模分析，提出了使用行生成算法將逆問題分為了主問題和子問題。通過對子問題中的整數(shù)規(guī)劃模型進行求解可以得到逆問題的最優(yōu)解。分析表明對于行生成算法中的子問題采用線性松弛法和本地搜索方法可以分別得到逆問題的上下界。數(shù)值實驗結(jié)果表明了線性松弛法得到的上界與最優(yōu)解較為接近，但求解效率并未有提升；啟發(fā)式算法在下界的獲取上則有著較為優(yōu)異的表現(xiàn)，同時求解速度較快，說明了該啟發(fā)式算法的有效性。

實際上，由于子問題是無容量限制設(shè)施選址的原問題，可以將求解該問題更為高效的近似算法用于其中，從而提高下界解的質(zhì)量和求解效率。子問題是最具有優(yōu)化空間的地方，可以結(jié)合啟發(fā)式算法一次找出多個目標(biāo)函數(shù)小于0的解，并一次性全部加入主問題中進行求解，從而提高求解效率。