几种基于匈牙利算法求解二次分配问题的方法及其分析比较

张惠珍;马良

【摘 要】二次分配问题(Quadratic assignment problem,QAP)属于NP-hard组合优化难题.二次分配问题的线性化模型和下界计算方法,是求解二次分配问题的重要途径.本文以二次分配问题的线性化模型为基础,根据现有QAP对偶上升下界计算方法中的具体操作,提出几种可行的QAP对偶上升计算新方法.最后,通过求解QAPLIB中的部分实例,深入分析其运行结果,详细讨论了基于匈牙利算法求解二次分配问题的对偶方法中哪些操作可较大程度地提高目标函数最优解的下界增长速度,这为基于匈牙利算法求解二次分配问题的方法的改进奠定了基础. 【期刊名称】《运筹与管理》 【年(卷),期】2010(019)001 【总页数】8页(P92-99)

【关键词】二次分配问题;下界;线性化;匈牙利算法 【作 者】张惠珍;马良

【作者单位】上海理工大学管理学院,上海,200093;上海理工大学管理学院,上海,200093 【正文语种】中 文 【中图分类】O224

1957年,KoopMans和BeckMann[1]在考虑经济活动选址问题时,不仅考虑活动位于特定位置的所需花费,而且考虑了活动之间的交互作用,首次提出二次分配问题

(Quadratic assignMentp rob lem,QAP)。该问题一般可描述为:给定n个设施和n个地点,三个n×n矩阵,即F=(fij)∈Rn×n,D=(dij)∈Rn×n,C=(cij)∈Rn×n,其中,fij表示设施i和j之间的流量,dij表示地点i和j之间的距离,cij表示设施i位于位置j的所需花费。要求给每个设施分配到一个地点,并使设备之间的总费用最小 QAP不仅综合了一大类组合优化问题的典型特征,而且被广泛应用于许多实际问题,如打字机键盘和控制面板设计,经济活动,调度问题,统计分析,及化学反应分析和考古数据排序等[2～4]。QAP问题通常在多项式时间内无法求解,属于NP-难解问题[5]。其主要原因是所谓的“组合爆炸”现象,求解时间随问题规模呈指数增长。一般而言,当问题规模时,很难在合理的计算时间内找到其最优解。

二次分配问题目标函数(1)中的二次项在一定程度上增加了问题的求解复杂度,通过一定方法将其二次项线性化,得到等价的二次分配问题线性化模型,不仅会使问题的求解复杂度得到一定降低,节省一定计算资源,而且能为原问题求得较优解或较优下界。半个多世纪以来,许多二次分配问题线性化方法已被提出,如:KaufMan和B roeckx线性化模型[6],夏勇和袁亚湘等在文献[7]和[8]提出的QAP线性化模型,Law ler线性化模型[9],Frieze2Yadegar线性化模型[10],A daMs2Johnson线性化模型[11]和Padberg2R ijal线性化模型[3]等。

本文在深入讨论二次分配问题线性化模型结构特征的基础上,探索了几种基于匈牙利算法的二次分配问题对偶上升求解法,通过求解二次分配基准问题库(QAPL IB[12])中的部分实例,对比分析其运行结果,说明了这些方法中哪些技术可取,哪些技术不可取。这为以匈牙利算法为基础的二次分配问题求解方法的研究奠定了基础,指明了一定方向。

通过引入新的变量yijkl=xijxkl(1≤i,j,k,l≤n),将QAP模型(1)中的二次项线性化,可得到如下与(1) -(2)等价的QAP线性化模型(简记为LQAP)

我们依照文献[13],将必然为“0”的变量yijil(1≤i,j,l≤n且j≠l)和yijkj(1≤i,j,k≤n且

i≠k)记为“无效变量(disallowed variables)”;“无效变量”以外的其它变量统称为“有效变量(allowed variab les)”。线性变量xij(1≤i,j≤n)记为“主导变量(leader variab les)”;满足yijkl=yklij的变量yijkl和yklij(1≤i,j,k,l≤n且i≠k,j≠l)记为“互补变量(comp leMentary variables)”。

由yijkl=xijxkl可知,LQAP中的四维数组可看作是二维排列矩阵x的克罗内克积(Kronecker Product),即y=xÆx。因而可将y看作是的矩阵n2×n2,并可将其分块为n2个的n×n分块矩阵。

其中,pij是变量xij的系数,pij=qijij+cij。我们将相应于“主导变量”、“有效变量”的系数分别称为“主导元素”和“有效元素”。

将y及Q分别用n-1条纵线和n-1条横线分成上述个子块,即y(i,j)和Q(i,j),1≤i,j≤n。不难看出,y和Q具有如下性质:

a)y和Q的分块矩阵中,每个子块有且仅有一个主导变量(元素),且各子块在分块矩阵y和Q中的相对位置与主导变量(元素)在该子块中的相对位置一致;

b)由y(i,j)=xijx及(4)可知,分块矩阵y每行每列有且仅有一个子块与解矩阵x相等,其余子块均为零矩阵。与解矩阵x相等的子块中,其主导变量必为“1”。这致使在保证子块q(i,j)中各元素非负的前提下,Q(i,j)中任一行(或列)的有效元素同时减(或加)一个常数,另一行(或列)的有效元素同时加(或减)该常数,不会改变原问题的解集及其目标函数值;

c)yijkl与其互补变量yklij(1l≤n且l)永远位于不同行,不同列的子块矩阵中。另外,在保证互补变量yijkl和yklij的系数非负的前提下,qijkl减去(或加上)一个常数,qklij加上(或减去)该常数,并不改变QAP的最优解;

d)将y看作是n2×n2的矩阵,yn2×n2的每行每列有且仅有一个元素取值为“1”,其余元素均为“0”。即yn2×n2中每行每列变量满足如下条件

这致使在保证Qn2×n2中各元素非负的前提下,Qn2×n2任一行(或列)的有效元素

加上(或减去)一个常数,并不改变原问题的最优解;

e)(4)是对yn2×n2中主导变量的约束,并可定义如下由主导变量构成的整数规划模型

匈牙利算法是基于线性分配问题最优解性质提出的,即如果将线性分配问题目标函数系数矩阵的每一行(列)的各元素都减去该行(列)的最小元素,得到一新的矩阵,则以新矩阵为系数矩阵的分配问题,其最优解与原问题最优解相同[14]。

根据QAP线性化模型的结构特征及匈牙利算法,可以提出多种求解QAP问题的对偶上升方法。他们在不改变原问题最优解的条件下,不断调整Q中各有效元素,通过求解L 2LAP和LAP(i,j)逐渐改善QAP最优解的下界值。经系列变化后,如果在分块矩阵Q中得到n个位于异行异列的子块,不仅在每个子块中存在n个位于异行异列且排列模式一致的零元素,而且各子块在Q中的相对位置与子块中主导元素在该子块中的相对位置一致(记具有这种零元素排列模式的矩阵Q为Q3),则已得最优解,解矩阵满足

2.1 Peter Hahn对偶上升下界求解法

Peter Hahn于1998年首次提出一种基于匈牙利算法的QAP对偶上升求解方法[13](简记为PHDF),该算法主要包括以下几个步骤 Step 1 求解LAP(i,j),1≤i,j≤n,并更新主导元素

除主导元素所在行列之外,对每个(n-1)×(n-1)的分块子阵Q(i,j),首先更新其有效元素的值:qijkl←qijkl+qklij,并将其所对应的互补元素置为“0”,即qklij←0。然后运用匈牙利算法求解LAP(i,j),并利用LAP(i,j)的目标函数值(记lij)更新主导元素pij:pij←pij+lij。

Step 2 求解L 2LAP,更新目标函数下界

利用匈牙利算法求解由主导变量构成的线性规划模型L 2LAP,将其目标函数值记为R,并更新原问题的目标函数值下界R′:R′←R′+R(R′的初值为0)。

记经过前述两步变化的目标函数系数矩阵为Q″。此时存在三种可能的情况:1)分块矩阵Q″中零元素的排列模式符合Q3中零元素排列模式,则已得最优解,R′为最优目标函数值;2)分块矩阵Q″中零元素的排列模式不符合Q3中零元素排列模式,R′的增量为0,则程序结束,R′为最优目标函数值的下界;3)分块矩阵Q″中零元素的排列模式不符合Q3中零元素排列模式,R′的增量不为0,则转Step 3。 Step 3 更新有效元素

在主导元素不为零的分块子阵Q(i,j)中,除主导元素所在行列外,利用pij/(n-1)更新其余有效元素,即qijkl←qijkl+pij/(n-1),并将其主导元素置“0”,pij←0。 Step 4 再次求解LAP(i,j)(1≤i,j≤n)并更新主导元素

该步与Step 1的不同在于求解LAP(i,j)(1≤i,j≤n)的先后顺序发生了变化。根据Step 2的计算结果,按照主导元素为“0”的子块优先计算,主导元素不为“0”的子块按从上到下,从左到右的顺序依次计算LAP(i,j)。然后,转Step 2。 2.2 新的对偶上升下界求解法

根据PeterHahn下界求解方法中的操作,我们提出如下问题:

(1)虽然在保证互补元素qijkl和qklij非负的前提下,相互转移他们之间的值,不会改变原QAP问题最优解,但上述Peter Hahn算法的每次循环都转移互补元素间的值,这对下界增长速度有一定影响吗?该操作对下界增长速度起到了积极的促进作用,或该操作在一定程度上对下界增长速度有一定的抑制? (2)除主导变量之外的有效变量也可构成如下线性分配模型

上述对偶上升下界求解过程中,某些步骤中是否可通过应用匈牙利算法求解L 2NLAP提高下界增长速度?

(3)由yijkl=yklij及yijkl=xijxkl可知,互补变量同时为“0”或同时为“1”。那么在对偶上升下界求解过程中,可通过qijkl←1/2(qijkl+qklij)和

qklij←1/2(qijkl+qklij)同时调整互补元素qijkl和qklij(1≤i补变量置于同等地位(经该操作后,qijkl和qklij对最优解的贡献量是相等的)。与PeterHahn对偶上升下界求解方法中调整互补元素的方法相比,哪种方法会较大程度地促进每次循环过程中的下界增量及下界逼近最优解的速度? 针对上述问题,我们提出几种新的对偶上升下界求解法,并将其分别命名为DF2I,DF2II,DF2III和DF2IV: DF2I对偶上升下界求解法:

DF2I与PeterHahn方法的不同在于其第一步,其余均相同。DF2I的第一步没有转移互补元素间的值,保留了每个中的原有数据结构特征。 DF2II对偶上升下界求解法:

DF2II与PeterHahn方法的不同在于其第一步和第三步,其余均相同。DF2II的第一步也没有转移互补元素间的值。DF2II的第三步更新有效元素时,除主导元素所在行列外,首先利用pij/(n-1)更新Q(i,j)中其余有效元素,即qijkl←qijkl+pij/(n-1),并将其主导元素置“0”,即pij←0;然后通过qijkl←1/2(qijkl+ qklij)和

qklij←1/2(qijkl+qklij)同时调整互补元素qijkl和qklij(1≤iDF2III是在Peter Hahn方法的基础上增加了求解L 2NLAP的操作。DF2III的前三步与Peter Hahn方法的相应操作相同,第五步与PeterHahn方法的第四步相同。DF-III的第四步是在PeterHahn方法的前三步结束后,首先通过

qijkl←1/2(qijkl+qklij)和qklij←1/2(qijkl+qklij)同时调整互补元素qijkl和qklij(1≤iDF2IV是在DF2II的基础上增加了求解L 2NLAP的操作。DF2IV的前三步与

DF2II的前三步相同,第五步与DF2II的第四步相同。DF2IV的第四步与DF2III的第四步相似,所不同的是,DF2IV的第四步求解L 2NLAP完毕后,应再次应用1/2(qijkl+qklij)同时调整互补元素qijkl和qklij。

本文算法用Matlab R2007a编制,实验环境为:3.40 GHz Pentium(R)4处理器,1GB内存,操作系统为W indow s XP Professional2002。我们选择QAPL IB中的部分实例来测试上述几种QAP对偶上升求解方法,有关结果见表1～表3。表1和表2给出37个的实例经PH 2DF,DF2I和DF2II求解的结果。如表1所示, Chr12a2Nug15,Rou152Nug16b,Nug172Nug18,Chr20a2Tai20b在由PH 2DF,DF2I和DF2II求解时所允许的最大循环次数分别为:100,50,30和20。 表3给出5个的实例经DF2III和DF2IV求解的结果。

表2所示37个求解算例中,仅有5个算例由PH 2DF求解的首次循环增量(FIncre)大于由DF2I和DF2II求解时的首次循环增量。这说明PH 2DF中第一步转移互补元素间的值,导致原问题Q(i,j)中的数据结构特征过早被破坏,对首次循环中的下界增量造成一定负面影响。由于在首次循环中,由DF2I和DF2II所计算的QAP下界增量明显大于由PH 2DF所计算的下界增量,这导致最终由PH 2DF所计算的绝大多数QAP实例的下界值小于由DF2I和DF2II所计算的下界值(如表1所示)。 由于DF2I和DF2II仅是首次循环操作相同,他们的首次循环增量相同。对比表2中由DF2I和DF2II所计算的除首次循环增量外的其余循环增量之和(R Incre)可知,DF2II中互补元素间值的转移方法更有助于QAP下界的增长。这导致经相同的循环次数,对大多数计算实例而言,DF2II的计算结果优于DF2I的计算结果(如表1所示)。但另一方面,这也给DF2II带来一定负面影响,如表2所示,DF2II计算耗时较大是其一大缺点。

综合考虑计算耗时和下界质量,从表1和表2,我们不难看出:DF2I在一定程度上优于PH 2DF和DF2II。这主要体现在:对大多数计算实例而言,DF2I不仅节省运算时

间,而且在一定程度上保证了解的质量。

对比表3和表1中相应实例的计算结果可知,DF2III和DF2IV与PH 2DF,DF2I和DF2II相比没有优势。事实上,DF2III和DF2IV中每次循环经计算LAP(i,j)和L-LAP后,L 2NLAP对QAP的下界增量贡献已很小,而且由于L 2NLAP规模较大,计算比较耗时,导致DF2III和DF2IV求解QAP的效率较低。

QAP有着广泛的应用背景,半个多世纪以来无数优秀的数学家为其倾注了大量精力,并已取得卓越的成果。尽管现已证明QAP是个NP—完备问题,是否存在多项式阶数算法(或有效算法)尚不可知,但我们相信,随着人们对优化理论和技术的深入研究,及现代计算机技术的快速发展,不仅会使既有的QAP求解方法不断改进,而且会出现一些既有理论依据又有实用价值的新算法。

【相关文献】

[1]KoopMans TC,BeckMannMJ.A ssignMen tp rob leMs and the location of econoMic activities[J].EconoMetrica,1957,25 (1):53276.

[2]Loio la EM,Abreu N MM,Boaventura2Netto PO,et al.A survey for the quad ratic assignMent p roblem[J].European JournalofOperationalResearch,2007,176(2):6572690. [3]PadbergMW,R ijalMP.Location,Scheduling,Design and Integer ProgramMing[M].Boston,K luwerA cadeMic Pub lish2 ers,1996.

[4]Pardalos PM,Rend l F,W o lkow iczH.The quad ratic assignMentp rob lem:a survey and recentdevelopMents[A].Pardalos PM,W o lkow iczH.Quadratic A ssignMentand Related ProbleMs[C].Providence,R.I:D IMACSSeries in D iscreteMathe2 Matics and TheoreticalComputer Science,AMericanMatheMatical Society,1994,16:1242. [5]Sahni S,Gonzalez T.P2comp lete app roxiMation p rob leMs[J].Journalof the A ssociation for CoMputingMachinery,1976, 23(3):5552565.

[6]KaufMan L,B roeckx F.An algorithMfor the quad ratic assignMent p rob leMusing Benders’decomposition[J].European JournalofOperationalResearch,1978,2(3):2072211. [7]X ia Y.IMp roved gilmo re2law ler bound for quad ratic assignMentp robleMs[J].Chinese Journalof EngineeringMatheMatics, 2007,24(3):4012413.

[8]X ia Y and Yuan Y.A new linearizationMethod for the quad ratic assignMentp roblem[J].Op tiMizationMethods and Soft2 ware,2006,21(5):8032816.

[9]Law ler E L.The quadratic assignMentp rob lem[J].ManageMent Science,1963,9(4):5862599.

[10]Frieze A M,Yadegar J.On the quadratic assignMentp rob lem[J].D iscrete App liedMatheMatics,1983,5(1):298.

[11]A daMsW P,Johnson TA.Imp roved linear p rogramMing2based low er bounds for the quadratic assignMentp rob lem[A]. Pardalos PM,W o lkow icz H.Quad ratic A ssignMentand Related Prob leMs[C].Providence,R.I:D IMACSSerieson D is2 creteMatheMatics and Theo reticalCoMpu ter Science,AMericanMatheMatical Society,1994,16:43275.

[12]Bu rkard R E,Karisch SE,Rend l F.QAPL IB 2a quadratic assignMen tp rob leMlib rary[J].Jou rnalof GlobalOp tiMization, 1997,10(4):3912403.

[13]Hahn P,Grant T.Low er bounds for the quadratic assignMentp rob leMbased upon a dual forMu lation[J].OperationsRe2 search,1998,46(6):9122922. [14]马良.基础运筹学教程[M].北京:高等教育出版社,2006.