基于DantzigWolfe分解和列生成算法的自存储收益管理大规模优化模型

基于Dantzig-Wolfe分解和列生成算法的自存储收益管理大规模优化模型摘要：自存储行业正处于高速发展阶段，每年10%的美国家庭租有一个自助存化Dantzig-Wolfe分解技术、分支定界技术提出了针对该模型的列生成算法，并用MATLAB程序实现了此算法。求解和比较表明，本文提出的列生成算法能准确解突出。关键词：自存储、收益管理、整数规划、列生成算法、D-W分解AbstractSelf-storageisaboomingindustryhasreceivedlittleacademicattentionfar.Inabout1outof10householdscurrentlyrentaself-storageunit.Thisconsidersaself-storagecompany,facingasetofreservationsforstorageunitsahorizonwithrevenuerewards.Thecompanyhastostoragerequeststobeacceptedtomaximizerevenue.Wemodelproblemasabinaryintegerprogramming,anditbyDantzig-Wolfedecompositionandcolumngeneration.WeimplementthisalgorithmbyTheexperimentsthatthisalgorithmcanefficientlyproblem,andcanmuchcomparedwithimprovedsimplexalgorithmanddualsimplexalgorithm.Keywords:Self-storageindustry,revenuemanagement,binaryintegerprogramming,columngenerationalgorithm,Dantzig-Wolfedecomposition1.引言20世纪60年代中期美国德克萨斯州，随后迅速蔓延至全美及欧洲。至2008年，美国已有30604家自助存储公51500大型自助存储仓库（self-storage国公民拥有6.89平方英尺自助存储仓库面积，将近10%美国家庭拥有一个租赁1995年为65%的增长[1]，也意味着自助存储已被越来越多的家庭认识、接受、并使用。上世纪80年代末，自助存储行业在英国出现，随后开始了其在欧洲的迅速2007869个。欧洲自助存储协会联盟（FederationofEuropeanSelfStorageAssociation,FEDESSA）和英国自助存储协会（SelfStorageAssociationoftheUnitedKingdom）的官方网站上都提供了该行业在欧洲的发展概况。不仅在Yellow，正在迪拜兴建一个全世界最大的自助存储仓库[1]。从以上数据及描述可以看出，自助存储行业虽然只出现了40多年时间，可是其市低迷及金融危机更赋予此行业发展的一个有利时机。本文研究的问题来源于Shurgard，好地控制了成本，因此，收入是该公司的主要关注点。Shuragrd公司的仓库都是收据公司的经验，大约90%的未来三个月内仓库预订需求可以提前知晓。有些客户大公司目前空置的仓库类别及各类别仓库数目限制。有收入。在现实生活中，该问题有五种情况：（1）只有一种仓库类别。每份订单只预订该仓库类别中的一个存储单元。这是最为简单的一种情况。（2）只有一种仓库类别。每份订单可能预订不只一个该仓库类别的存储单元。这种情况在一些中小型自助存储公司中很典型。（3）有不只一种仓库类别。每份订单总共只预订某种仓库类别的一个存储单元。这是最普遍的一种情况。在Shurgard，大部分家庭客户和一些商业客户的订单都属于这种情况。（4）有不只一种仓库类别。每份订单可能预订多个存储单元，并且这些存储单元同属于一类仓库类别。许多商业客户的订单属于这种情况。（5）有不只一种仓库类别。每份订单可能预订多种仓库类别中的多个存储单元。这种情况中订单预订没有任何限制，是最一般的情况。由于第3种情况和第4种情况中的订单只能预订一类仓库，因此，可以将情况3分解为多个属于情况14分解为多个属于情况21和情况23和情况4所属的问题即可相1和情况2所属的问题可用多项式算法解决，将其转化为最小成本流问题，算法复杂度为O((nn)2)，从而解决了问题1和问题也解决了问题3和问题4。由于这四种情况并非本文研究问题，这里对其不再赘述。本文主要研究情况55，算法复杂度为O(nm+1)表示所有类别仓库数量总和。可以看出，当仓库数量较多的时候，该问题的动态规划算法将会耗费大量的求解时间[2]。本文提出的列生成算法将从另外一个角度对问题5进行建模并对其作更有效地求解。2.模型2.1模型符号定义在建立问题的模型之前我们首先用数学标号对问题进行数学定义问题的一些符号如下表所示：t时t=间单位。j订单编号。j=2,…,ni仓库类别总数。i=2,…,KM所有仓库总数。第i类仓库数量。第j份订单预订的第i类仓库的数量。Sj第j份订单货物的开始储存时间。Cj第j份订单货物的结束储存时间。每个第i类仓库每时间单位的收费价格。vj接受第j份订单产生的收益。值得注意的有以下几点：（1）（2）

K∑i=mi=1Kj≤因而∑j≤,i=1

j=2,…,n（3）Sj和Cj分别表示第j份订单的储存的开始时间和结束时间，因此，对于订单1要求其货物于第2个时间单位到第4个时间单位之间被储存，则S1=2,C1=4,因此，其被储存时间为（4-2+1）=3个时间单位。K K（4）

vj=∑iSj⋅(Cj−Sj+),i=1

j=1,2,…,n.令Wj=∑i=1vj=Wj⋅(Cj−Sj+2.2原模型及D-W分解模型(1)问题的整数规划模型及松弛模型为了建立此问题的整数规划模型，本文假设了两类变量：第一类变量：xj，其为0-1变量。若xj=1，则表示订单j最终被接受；否则xj=0。第二类变量：xjt，其也为0-1变量。对于任意

j2,…,n}，t∈{Sj,Sj+,Sj+2,…,Cj}

。若xjt=1，则表示t时刻j订单

t∈{Sj,Sj+Sj+2,…,Cj}=0。显然，当xj=1时，对任意xjtxj=0时此，

Cj∑xjt=(Cj−Sj+⋅xjt=Sj则问题的整数规划模型（IP）如下：n n CjZIP=max

∑j=1

vjxj

=max

∑∑j=1t=Sj

WjxjtCj∑xjt−(Cj−Sj+⋅xj=0,t=Sjn

j=2,…,n

（1）∑xjt≤,j=1

i=K;

t=T

（2）

（3）将约束（3）松弛到[0,1]，则原问题可松弛为如下线性规划问题(LP)：n n CjZLP=max

∑

vjxj

=max

∑∑j=1t=Sj

WjxjtCj∑xjt−(Cj−Sj+⋅xj=0,t=Sjn

j=2,…,n

（4）∑xjt≤,j=1

i=K;

t=T

（5）0≤≤1,0≤≤1

（6）其中约束（4）表示若订单j被接受，则其对应的所有变量均为1；否则，其对应的所有变量均为0。约束（5）表示仓库资源约束，即任意一时刻任意一类仓库的使用数不能超过这一类仓库的数量。(2)D-W分解原理及问题的D-W分解模型令约束（5）的系数矩阵为S，X为xjt组成的向量，注意到，S为块对角矩阵，因此，约束（5）可拆分为如下T组约束条件：⎧Xj11≤m1≤m,Xj2 i2⎨2⎪SSX ≤m,i

titi⋮t=T,i⎩T其中，Xjt,t2,…,T}表示所有满足条件：Sj≤t≤Cj的xjt组成的向tt时刻可以被加工的订单所组成的变量集 Sij 合。表示对应变量的系t数。另外，原模型的约束条件个数为n+K*T个，因此，若订单数和时间范围量的时间。结合以上两点，本文采用分解（Dantzig-Wolfe分解）对原模型进行分解变形，使之相对于原模型有以下几点改进：（1）变形后的分解模型将会使约束条件大大减少，为n+T个，虽然分解后变量个数将会比原模型多出许多，然而，本文采用的列生成算法可以很好地解决变量数目巨大这一问题；（2）由于可将原问题中的系数矩阵S分解为T组约束条件（每组由K分解后，列生成算法的子问题为T个独立的具有K个约束条件的线性规划问题，这使得子问题的求解难度大大降低。（3）经过分解后的模型下界与原模型相比要更为收紧，因此，只需进行有限步迭代即可获得原问题松弛模型的最优解，节省了求解松弛模型的时间。分解原理的描述如下：考虑任意LP问题，称为原问题：x∈S

(7)其中A为m*n矩阵假设S是有界多面体令1,…,xP是S的极点集合则有：∑p ∑p p∈

λxp,

λλ

≥

（8）problemorMP）：

P∑

p(cxppP∑pPpp)λ=b

（9）∑λp=1λp≥n问题（7）和（9）是等价的，它们有同样的最优解。问题(4.9)的任意可行解λ相应问题（7）的可行解x,反之亦然。这就是Dantzig-Wolfe分解原理，简称分解。在分解中，用单一加权约束代替了多面体集约束减少了约束数目，但因这个多面体的顶点可能非常多（Cm种多列的LP分解的优势体现在它收紧了原问题的下界，而且只需有限步迭代即可获得原问题的最优解。n根据以上对D-W分解原理的描述，约束（5.1）至（5.T）中的每组约束与约束（6）的可行域分别构成一个多面体，设为,P2,…,PT。它们分别代表着tt t时刻仓库资源约束的可行域多面体。令Xl(l2,…,L})表示多面体Pt的第t tl lt个顶点，其中表示相对应于t时刻的多面体的顶点数。则Pl lt

内的点Xt可表示Xt=

Ltl l ll=1

λtXt,其中

Lt∑l=1

λt=1，λt≥0因此，问题LP可被转化为以下模型，称为DWLP：T ⎛n ⎞llZDWLP=maxll

∑∑⎜∑

Wjxjt⎟

λtT

t

l⎝

j⎠∑∑jt t j j j

（10）tl

λl−(C−S

+⋅x

=0,

j=2,…,nl∑l=

t=T

（11）0≤λl≤1,0≤

x≤1

（12）t j至此，原问题被转化为分解后的线性规划问题。可以看到，原问题有tn+K*T个约束条件，而经过分解后的问题约束条件为n+T中变量为λl），t不过用列生成算法来求解这个问题，变量的个数将不再成为困扰。3.算法基具化可以为求解本文的问题模型提供更多的启发。(1)列生成算法基本思想与原理列生成算法（ColumnGenerationAlgorithm）的基本出发点是研究这样一类大规模线性规划问题：x≥0

（13）nx,A,b是有合适维数的矩阵，在这些问题里变量数n远远超过约束数m，以至于计算机不能对约束矩阵进行有效的存储和计算。这种例子如分解的主问题里，每一列对应一个基可行解，即可行域的顶点，因为有Cm个顶点，当n很大时，即使m较小，总的极点数即主问题列数（变量）也是一个非常巨大的数，处理这类问题，可以使用列生成算法。它只是用主问题列的一些子集所构成的问题，称为限制主问题（Restrictedmasterproblem,RMP）去求解原问题，这就是列生成算法的基本思想。事实上，m，这样就可以使用较少的列获得原问题的最优解，而其它的列只有当需要的时候才去产生。n列生成算法的出发点是这样的，由于所研究问题的大规模特性，使得不可法获得这种大规模问题的解。列生成算法框架如图（4-1）所示。假设在求解过程的某一步，一个限制主问题（RMP）A'x=b，它是原问题（主问题）的一个缩小问题，仅包含原问题系数矩阵A的部分列。通过解这个限制主问题，获得单纯形乘子（simplexπ=(π1,…,πm)。然后利用限制主问题的解信数：moj=cj−∑πiaij≥0,i=1

,n则满足对偶可行性，表示已获得原问题的最优解。但假使存在某一列s，使σj进行直到获得原问题的最优解。可是问题在于原问题的列非常多，不可能用枚举方法找到这样的列，因此需构造一个子问题（subproblem,（pricingproblem）来获得这样的列或判断原问题解的最优性。典型地，通过使用这种交互方法，只需有限步迭代即可获得原问题的最优解。综上，列生成算法一般步聚如下：1：用启发式算法（或其它方法）生成初始解，每一个子问题至少生成2：解限制主问题，求得单纯形乘子π。3：解相应价格问题（对偶问题），判断最优性，如果得到最优解，终止；否则，转4。4：选择具有最小负判别数，生成添加列。5：返回2。获得原问题最优解获得原问题最优解解限制主问题：min∑c'jxjst. A'x=bx≥0'为单纯形乘子将列s加到限制主问题：A'=A+s解子问题：m mβ=min(c'j−∑πiaij)=c's−∑πiaisi=1 i=1β≥0?NY图1：列生成算法框架(2)列生成算法在块对角状大规模线性规划上的应用列生成算法最早是被用来求解块对角状大规模线性规划问题（DantzigandWolfe,1961;Gilmoreand式如下：X1+C2X2+⋯+CpXp+X2+⋯+Xp=

（14）

⋮

Xp

===Xi≥i=p 10由分解原理，上述LP问题对应的主问题（MP）为：p ∑∑ i i ijz=Xjp j∑∑(=j

（15）

j

≥

j则提出问题列生成算法步骤如下：i i1：利用两阶段算法，对每个i=1,…,p求得两个初始极点X2i iπs=|π0

]到初始的RMP。解得：1，1，λ

p1，λp2

及单纯形乘子πs

，分解πs为两部分：πs=|π0]，其中π对应组合约束，π0对应p个凸组合约束，k=k0

+1，这里k0=2是初始极点数。2：解子问题：=−π)Xi−π0iXi=Xi≥0对所有i=p，获得极点

Xk，其中πi是π

的第i个分量。如果i 0 0ui≥0,i=1,…,p，那么RMP得最优解，也就是原问题的最优解，算法停止；否则，令us=:ui<0}，转3。3：产生添加列：s s⎛AXks s=⎜ ⎟⎝ ⎠其中，es是p维单位向量，并且它的s分量为1。4：解新的RMP，得,…,λ1,…,λ,

及π'，返回2。p plp s(3)自助存储优化问题的列生成算法根据以上两部分对列生成算法的描述，回到自助存储优化问题的模型，对DWLP中的T个多面体，首先对每个多面体取两个顶点，X1,X2,X1,X2,…,X1,Xλ1,λ2,λ1,λ2,…,λ1,λ1 1 2 2 T T

1 1 2 2 T TDWLP构成初始限制主问题（n+T个约束条件，因此，用单纯形法即能求解。得到单纯形乘子π=j|πt]，其中，πj表示对应约束（4.10）第j个约束条件的单纯形乘子；πt表示对应约束（4.11）第t个约束条件的单纯形乘子。由于原问题的资源约束矩阵可以分解成T个独立的T个独立的子问题。由于每个子问题是一个具有K个约束条件的线性规划，而解T个具有K个约束条件的线性规划的难易程度及耗费时间要大大小于解一个具有SP(t)如下：(SP(t)):n nxσt=Wjxt−∑πjxt−πtn

j=1. ∑Sjxt≤i,j=1x

i=1,2,…,Kt为常数，因此实际上子问题是求目标函数右边部份的前两个表达式之差。由于子问题只有K个约束条件，而在现实生活中，K往往是一个不大的数。这T个子问题可以很容易地由单纯形法求解。取στ

=t}，若στ≤0，则所求的解即为DWLP的最优解；否则，生成相应的列加入RMP中继续求解，直到满足以上的条件，则所求解为最优解。由于DWLP的最优解不一定是IP问题的最优解，甚至不是IP问题的可行解，因此，接下来需要对非整数变量进行分支定界。本文采用的分支策略如下：ss j对求得的DP的最优解*的第一类变量进行判断取其中一变量x*使其满足以下条件：|x*−05=max(|x*−05|),j=,…,n，对其进行分支定界，之所以这样取的原因是因为此变量最不“整数性（罗守成20。将s分为两支s=；s=，分别加进原问题约束条件中。事实上，当s=0时，与s相对应的第二类变量xt均为0；当s=1时，与s相对应的第二类变量均为ss j由以上分析，本文得出所研究问题的列生成算法步骤：Stepvj按降序排列，利用贪婪算法求得一原问题的可行解，作为原问题的初始下界。Step2：对每个多面体各取两个顶点X1,X2,X1,X2,…,X1,X1 1 2 2 T T由于变量均小于等于1，且Sij≤，因此，初始两个顶点可以选取如下：1X1=11X2={0,1,0,…,0}1⋮TX1=TTX2={0,1,0,…,0}T以上向量的维数与顶点向量维数相同（事实上，每个多面体各取一个初始顶将顶点及相对应的λ变量代入DWLP形成初始限制主问题如下：T 2⎛n ⎞llZRMP=maxll

∑∑⎜∑

Wjxjt⎟

⋅λttT 2

l⎝

j⎠∑∑jt t j j jtll2l

λl−(C−S

+⋅x

=0,

j=2,…,n∑l=1=

t=T0≤λl≤1,0≤

x≤1t j求得初始解及单纯形乘子π=j|πt]。Step3：令πj表示对应约束（4.10）第j个约束条件的单纯形乘子；πt表示对应约束（4.11）第t个约束条件的单纯形乘子。同时求解T个子问题，得到στ=t}。其对应的子问题的解为SPXτ。Step4：判断解的最优性。若στ≤0，则所求的解即为DWLP的最优解；否τ则，生成对应的列如下加入RMP中重新求解，其对应的变量为λ3：ττ(0,…,0,τ

τ,0,…,0,eτ

)T，其中，e

是一个T单位向量，其第τ个元素为1；之前的部分是一个n维向量，SPXτ对应于τ时刻的订单所在的行。jStepRMP，直至求得DWLP的最优解λ*和X*j

（第一类变量j最优值）。若X*为整数，则其亦为原IP问题第一类变量的最优解，将λ*还原成jtX*，原问题解毕，程序结束；否则，转6tStep6：对

X*进行判断，取其中一变量

x*，使其满足以下条件：js|x*−0.5max(|x*−0.5|),j=2,…,n，对其进行分支定界：x*=0；x*=1。jss j s sStep7：将两分支分别加入原问题IP，转2。4.应用规模较小的模型为例子，并详细说明如何使用上节所述算法求解该模型。此模型背景如下：有一自助存储公司，其总共有2类仓库，每类仓库数量分别有2个和35个订单需要公司来决定作决策，这5个订单所包含的时间段为第1个时间单位至第5个时间单位。因此，在此模型中，T==

K===3

，并且本文假设==5，Sj,Cj,Sij都根据相应的约束在excel里随机生成，如下表：订单jSjCjjS2jWjvj1231224510483450154141251202由上表可以得到此问题整数规划模型如下：Z

=+8x2+10x3+56x4+20x5=14(x12+x13)+24+x25)+34+x35)41+x42+x43+x44)51+x52)x1213−1=0x2425−2=0x34+x35−2x3=0x41424344−4=0x5152−5=0x41≤2x41+51≤3

tx1242≤2x12+2x42+2x52≤3x1343≤13+43≤3x2444≤x34+44≤3x25≤x35≤3

t=t=t=t=xj

j2,…,5},t对上述模型进行松弛后，对每个多面体各取初始顶点，为了计算简便，本例对每个多面体各取一个初始顶点：

X1=(x,x

)T=，1 41 51X1=(x,x,x=，

X1=(x,x=，2

3 X1=(x

,x,x=，=(x

,x分别对应着λ1,λ1,λ1,λ1,λ14

5

1 2 3 4 5五个变量，代入以上模型松弛后可得该问题的初始RMP为：1 2 3 4 5maxZRMP=14(λ1+λ1+λ1)+1 2 3 4 52 3 1λ1+λ1−2x2 3 14 5 2λ1+λ1−2x4 5 2t jλ1=1,0≤t j

−2=01 4λ1−2x=1 4−2=0≤2,…,5},j用单纯形法求解，得最优值：Z*RMP=50,单纯形乘子π=(0,0,…,0)。分别解5个子问题：SP1=max14x

+10x

SP2=max14x

+14x

+10x41 51s.t.

s.t.

12 42 52≤2+≤22x41+2≤3

2+2x42+2≤30≤x41,≤1

0≤,,≤1SP3=max14(x

+x)

SP4=max4x

+5x

+14x13 43s.t.

s.t.

24 34 44+≤2+≤22(+)≤3+2≤30≤,≤1

0≤x24,x34,x44≤1SP5=max4xs.t.

25+5x35≤2≤30≤x25,x35≤1解得=SP2*=21,SP3*=21,SP4*=23,SP5*=9,其中最大值为SP4*=2