第5章整数规划

整数规划解决决策变量是整数值的线性规划和非线性规划问题，是在有限个可供选择的方案中寻找满足一定标准的最好方案。整数规划自1958年由R.E.戈莫里提出割平面法之后形成独立分支，60多年来发展出很多方法。目前比较成功又流行的方法是分枝定界法和割平面法。01规划在整数规划中有重要的地位，有界变量的整数规划与01规划等价，同时许多实际问题（如背包问题、送货问题、指派问题等）都可以归为这类问题。整数规划的应用范围也极其广泛。它不仅在工业和工程设计、科学研究方面有许多应用，在计算机设计、系统可靠性、编码和经济分析等方面也有新的应用。
5.1本章内容
本章主要介绍以下内容。
（1） 整数规划的基本概念。
（2） 实际例子的建模方法。
（3） 割平面法及其MATLAB实现。
（4） 分支定界法及其MATLAB实现。
（5） 01整数规划及其MATLAB实现。
解整数规划最典型的做法是逐步生成一个相关的问题，称为原问题的衍生问题。对每个衍生问题又伴随一个更易求解的松弛问题（衍生问题称为松弛问题的源问题）。通过松弛问题的解来确定其源问题的归宿，即确定源问题应被舍弃还是再生成一个或多个它本身的衍生问题来替代它。随后再选择一个尚未被舍弃或替代的原问题的衍生问题。重复以上步骤，至不再有未解决的衍生问题为止。
整数线性规划问题与线性规划问题相类似，只是一些变量要求为整数。要求部分或全部决策变量必须取整数值的规划问题称为整数规划。不考虑整数条件，由余下的目标函数和约束条件构成的规划问题称为该整数规划问题的松弛问题。若松弛问题是一个线性规划，则称该整数规划为整数线性规划。如果没有连续的变量，则是整数规划问题。
例如，给定矩阵A、B和向量b、c、d，则有



min c′x+d′y




使得



Ax+By=b

x,y≥0

x是整数







这是一个混合整数线性规划问题，即不是其中所有变量均需要取整数值。
给定如上整数线性规划问题，对应的松弛问题即为



min c′x+d′y




使得



Ax+By=b


x,y≥0




整数线性规划数学模型的一般形式为



max（或min）z=∑nj=1cjxj




使得



∑nj=1aijxj≤(或=,或≥)bi,i=1,2,…,m

xj≥0,j=1,2,…,n

x1,x2,…,xn是整数




整数线性规划问题可以分为下列几种类型。
（1） 纯整数线性规划（pure integer linear programming）。指全部决策变量都必须取整数值的整数线性规划，有时也称为全整数规划。
（2） 混合整数线性规划（mixed integer linear programming）。指决策变量中有一部分必须取整数值，另一部分可以不取整数值的整数线性规划。
（3） 01型整数线性规划（zeroone integer linear programming）。指决策变量只能取值0或1的整数线性规划。
5.2建模方法
本节将介绍用于解决线性规划问题的一些建模方法。因为没有一个系统的方法来公式化地离散最优化问题，所以设计一个好的模型显得尤为重要。
5.2.1二元选择
可以运用一个二进制变量x，根据两个选项的选择将x置为0或1。
若变量只能取值0或1，称其为01变量。01变量作为逻辑变量（logical variable），常被用于表示系统是否处于某个特定状态，或者决策时是否取某个特定方案。例如



x=1,当决策取方案P时


0,当决策不取方案P时（即取时）




当问题含有多项要素，且每项要素均有两种选择时，可用一组01变量来描述。一般地，设问题有有限项要素E1,E2,…,En，其中每项Ej有两种选择Aj和j(j=1,2,…,n)，则可令



xj=1,若Ej选择Aj


0,若Ej选择jj=1,2,…,n




在应用中，有时会遇到变量可以取多个整数值的问题。这时，利用01变量是二进制变量（binary variable）的性质，可以用一组01变量来取代该变量。例如，变量x可取0与9之间的任意整数时，可令



x=20x0+21x1+22x2+23x3≤9




其中，x0,x1,x2,x3皆为01变量。
01变量不仅广泛应用于科学技术问题，在经济管理问题中也有十分重要的应用。
例5101背包问题。
有n个物品，第j个物品的重量为wj，价值为cj。一个背包所能承受重量的上限为K。如何选择物品使得总价值最大？这里假定不能选取部分物品，这种情况称为01背包问题（示意图如图51所示）。


图5101背包问题示意图


解： 为了给这个问题建模，可以引入二进制变量xj。如果物品j被选择，那么xj=1，否则xj=0。这个问题可以公式化为



max∑nj=1cjxj




使得




∑nj=1wjxj≤K

xj∈0,1,j=1,2,…,n




例52固定费用问题。
有3种资源(A,B,C)被用于生产3种产品，资源量、产品单件可变费用及售价、资源单耗量及组织3种产品生产的固定费用见表51。要求制订一个生产计划，使总收益最大。


表51资源及产品的相关数据



项目单耗量

产品Ⅰ产品Ⅱ产品Ⅲ资源量

资源A248500
资源B234300
资源C123100
单件可变费用456
固定费用100150200
单件售价81012


解： 总收益等于销售收入减去生产上述产品固定费用和可变费用之和。建模的困难主要在于事先不能确切知道某种产品是否生产，进而不能确定相应的固定费用是否发生。下面借助01变量解决这个困难。
设xj是第j种产品的产量，j=1,2,3； 再设



yi=1,若生产第j种产品（即xj>0）


0,若不生产第j种产品（即xj=0）j=1,2,3




则问题的整数规划模型是



max z=(8-4)x1+(10-5)x2+(12-6)x3-100y1-150y2-200y3




使得


2x1+4x2+8x3≤500

2x1+3x2+4x3≤300

x1+2x2+3x3≤100

x1≤M1y1

x2≤M2y2

x3≤M3y3

xj≥0且为整数,j=1,2,3

yj=0或1,j=1,2,3




其中，Mj为xj的某个上界。例如，根据第3个约束条件，可取M1=100，M2=50，M3=34。
如果生产第j种产品，则其产量xj>0。此时，由约束条件xj≤Mjyj，知yj=1，相应的固定费用在目标函数中将被考虑。如果不生产第j种产品，则其产量xj=0。此时，由约束条件xj≤Mjyj可知，yj可以是0，也可以是1。但yj=1不利于目标函数z的最大化，因而在问题的最优解中必然是yj=0，从而相应的固定费用在目标函数中将不被考虑。
5.2.2强制约束
离散最优化问题的一个共同特点是一些变量之间相互不独立。特别地，假设只能在选择了B的情况下才能选择A，可以引入二进制变量x(y)，如果选择A(B)，则x=0(y=0)。那么A与B之间的关联就可以表示为



x≤y




例53设施选址问题。
如图52所示，设有n个可能建设设施的地点，m个需要的用户。在j地建设设施需要花费cj，为用户i提供设施j需要花费dij。如何建设设施使得总花费最少？


图52设施选址问题示意图



解： 对每个地点j，设二进制变量yj，如果j处建设了设备，那么yj=1，否则yj=0。设二进制变量xij，如果给用户i配备设备j，那么xij=1，否则xij=0。设施选址问题公式化为（FL）



min∑nj=1cjyj+∑mi=1∑nj=1dijxij




使得




∑nj=1xij=1,i

xij≤yj,i,j

xij,yj∈0,1,i,j




如果地点j处没有设备，即yj=0，那么用户不可能配备j处的设备，则xij=0。
5.2.3变量之间的关系
约束条件



∑nj=1xj≤1




其中所有的变量都是二元的。意味着最多只有一个xj能为1。类似地，约束条件



∑nj=1xj=1




那么有且只有一个变量为1。
5.2.4析取约束
设x为一个非负的决策变量。假设有两个约束条件a′x≥b，c′x≥d，其中a′和c′均为非负。
设二进制变量y，有以下约束条件



a′x≥yb

c′x≥(1-y)d

y∈{0,1}


更一般地，假设有m个约束条件a′ix≥bi，i=1,2，…,m，对任意的i，有ai≥0，并且至少其中k个满足。引入m个二进制变量yi，i=1,2，…,m，有以下约束条件



a′x≥biyi，i=1,2，…,m

∑mi=1yi≥k

yi∈0,1，i=1,2，…,m




5.2.5值的约束范围
假设变量x的值在集合{a1,a2，…,am}中取得。可以引入m个二进制变量yj，j=1,2，…,m，有约束条件



x=∑mj=1ajyj

∑mj=1yj=1

yj∈{0,1}




5.2.6分段线性成本函数
成本函数为定义在［a1,ak］上的分段线性函数f(x)。区间［a1,ak］上的x可以表示为



x=∑ki=1λiai




其中，λ1,λ2,…,λk非负且和为1。
但这样x的表示不是唯一的，假定至多两个相邻的λi是非零的，这样任何的x∈［ai,ai+1］可以表示为x=λiai+λi+1ai+1且



f(x)=∑ki=1λif(ai)




引入二进制变量bi，i=1,2,…,k，公式化至多两个相邻的λi非零的条件，得到这个问题的公式化形式为



min∑ki=1λif(ai)




使得




∑ki=1λi=1

λ1≤y1

λi≤yi-1+yi,i=2,3,…,k-1

λk≤yk-1

∑k-1i=1yi=1

λi≥0

yi∈{0,1}





5.3整数规划的例子
例54某服务部门各时段（每2h为一时段）需要的最少服务员人数见表52。按规定，服务员连续工作8h（即4个时段）为一班。现要求安排服务员的工作时间，使服务部门服务员总数最少。


表52各时段需要的最少服务员人数




时段12345678
最少服务员人数10891113853

解： 设在第j时段开始时上班的服务员人数为xj。由于第j时段开始时上班的服务员将在第（j+3）时段结束时下班，故决策变量只需要考虑x1、x2、x3、x4、x5。
问题的数学模型为



min z=x1+x2+x3+x4+x5




使得



x1≥10
x1+x2≥8
x1+x2+x3≥9
x1+x2+x3+x4≥11
x2+x3+x4+x5≥13
x3+x4+x5≥8
x4+x5≥5
x5≥3
x1,x2,x3,x4,x5≥0,且均取整数值




这是一个纯整数规划问题。
例55现有资金总额为B。可供选择的投资项目有n个，项目j所需投资额和预期收益分别为aj和cj(j=1,2,…,n)。此外，由于种种原因，有3个附加条件： 第一，若选择项目1，则必须同时选择项目2，反之则不一定； 第二，项目3和项目4中至少选择一个； 第三，项目5、项目6和项目7中恰好选择两个。应当怎样选择投资项目，才能使总预期收益最大？
解： 每个投资项目都有被选择和不被选择两种可能，为此令



xj=1,对项目j投资


0,对项目j不投资j=1,2,…,n




这样，问题可表示为



max z=∑nj=1cjxj




使得



∑nj=1ajxj≤B
x2≥x1
x3+x4≥1
x5+x6+x7=2
xj=0或1，j=1,2,…,n




这是一个01规划问题。其中，中间3个约束条件分别对应3个附加条件。
例56工厂A1和A2生产某种物资。由于该种物资供不应求，故需要再建一家工厂。相应的建厂方案有A3和A4两个。这种物资的需求地有B1、B2、B3、B4 4个。各工厂年生产能力、各地年需求量、各工厂至各需求地的单位物资运费cij见表53。


表53各工厂及各需求地的相关数据



Ai
单位物资运费cij/（万元/千吨）

Bj
B1B2B3B4生产能力/

（千吨/年）


A12934400
A28357600
A37612200
A44525200
需求量/（千吨/年）350400300150
工厂A3和A4开工后，每年的生产费用估计分别为1200万元和1500万元。现要决定应该建设工厂A3还是A4，才能使今后每年的总费用（即全部物资运费和新工厂生产费用之和）最少。
解： 这是一个物资运输问题，其特点是事先不能确定应该建A3和A4中的哪一个，因而不知道新厂投产后的实际生产费用。为此，引入01变量



y=1，若建工厂A3


0，若建工厂A4




再设xij为由Ai运往Bj的物资数量(i,j=1,2,3,4)，单位是千吨； z表示总费用，单位是万元。
问题的数学模型为



min z=∑4i=1∑4j=1cijxij+［1200y+1500(1-y)］




使得


x11+x21+x31+x41=350
x12+x22+x32+x42=400
x13+x23+x33+x43=300
x14+x24+x34+x44=150
x11+x12+x13+x14=400
x21+x22+x23+x24=600
x31+x32+x33+x34=200y
x41+x42+x43+x44=200(1-y)
xij≥0，i,j=1,2,3,4
y=0或1（51）




上述数学模型中，目标函数由两部分组成，和式部分为由各工厂运往各需求地的物资总运费，加号后的中括号部分为建工厂A3或A4后相应的生产费用。约束条件式（51）中的前8个约束条件为供需平衡条件。第7个和第8个约束条件中含01变量y。若y=1，则表示建工厂A3。此时，第7个约束条件就是对工厂A3的运出量约束。再由第8个约束条件，必有x41=x42=x43=x44=0； 反之，若y=0，则表示建工厂A4。
显然，这是一个混合整数规划问题。
5.4问题的公式化
因为计算的复杂度是根据变量数n，限制条件数m呈指数增长的，要更好地公式化一个问题，就要尽可能减少使用的变量数和限制条件数。
这个问题就从研究整数线性规划问题的松弛问题入手。松弛问题的解也是整数线性规划问题的解。
在例53提到的设施选址问题中，可以考虑以下的等价的公式化： 



min∑nj=1cjyj+∑mi=1∑nj=1dijxij




使得



∑nj=1xij=1，i
∑nj=1xij≤myj，j
xij,yj∈0,1，i,j




这种公式化叫作aggregate facility location formulation（AFL）。
可以看到限制条件∑nj=1xij≤myj使得当yj=0时，xij=0，如果yj=1，那么xij可以等于1。也就是说这个限制条件与之前得出的xij≤yj,i,j是等价的。所以两组公式化均可解决设施选址问题。原来的公式化得到的结果有m+mn个限制条件，而新得到的公式只有m+n个限制条件。
为了比较两种公式化，考虑与它相对应的线性规划松弛问题。对于松弛问题，整数的限制条件xij,yj∈0,1改为0≤xij≤1，0≤yj≤1。
两个松弛问题的可行解分别为



PFL=(x,y)∑nj=1xij=1,i

xij≤yj,i,j

0≤xij≤1，0≤yj≤1

PAFL=(x,y)∑nj=1xij=1,i

∑nj=1xij≤myj,j

0≤xij≤1，0≤yj≤1




显然，PFLPAFL。也就是说，FL公式化方法对应的松弛问题的可行解比AFL方法的更加接近整数线性规划问题的解，但是AFL公式化方法限制条件的数量有非常明显的减少。
那么什么是整数线性规划问题理想的公式化呢？设T={x1,x2,…,xk}是某个整数规划问题的可行整数解的集合。假设T是有限的，考虑T的凸包



CH(T)=∑ki=1λixi∑ki=1λi=1,λi≥0,xi∈T




集合CH(T)是一个极值点为整数的多边形，且任何线性规划松弛问题的可行解集合P满足CH(T)P。如果我们准确地知道CH(T)，即可以将CH(T)表示为CH(T)={x|Dx≤d}，就可以通过寻找线性规划问题



min c′x




使得



x∈CH(T)




的极值解来解决线性整数线性规划问题



min c′x




使得



x∈T




理想的公式化得到的线性规划松弛问题的解就是整数线性规划整数可行解的凸包，但这是很难做到的。为得到较好的公式化，希望能够得到尽量接近CH(T)的多边形。
整数规划问题公式化的质量可以由松弛问题解与整数可行解的凸包CH(T)来决定。
一般情况下，松弛问题的最优解不会刚好满足变量的整数约束条件，因而不是整数规划的可行解，自然也不是整数规划的最优解。此时，若对松弛问题的这个最优解中不符合整数要求的分量简单取整，则所得到的解不一定是整数规划问题的最优解，甚至也不一定是整数规划问题的可行解。
例57考虑下面的整数规划问题



max z=x1+4x2




使得



-2x1+3x2≤3
x1+2x2≤8
x1,x2≥0




解： 图53中四边形OBPC及其内部为松弛问题的可行域，其中整数格点为整数规划问题的可行解。根据目标函数等值线的优化方向，从直观可知，P点（x1=18/7，x2=19/7）是其松弛问题的最优解，其目标函数值z=94/7。在P点附近对x1和x2简单取整，可得四点： A1、A2、A3和A4。其中，A1和A2为非可行解； A3和A4虽为整数可行解，但不是最优解。本例整数规划的最优解为A*点（x1=4，x2=2），其目标函数值z=12。
由于整数规划及其松弛问题之间的上述特殊关系，像例57中先求松弛问题最优解，再用简单取整的方法虽然直观简单，却并不是求解整数规划的有效方法。


图53四边形可行域示意图



5.5割平面法
考虑整数规划问题



max z=∑nj=1cjxj




使得


∑nj=1aijxj=bi,i=1,2,…,m

xj≥0,j=1,2,…,n

xj取整数,j=1,2,…,n（52）




设其中aij(i=1,2,…,m;j=1,2,…,n)和bi(i=1,2,…,m)皆为整数（若不为整数时，可乘上一个倍数化为整数）。
割平面法的基本思路如下。
（1） 解对应的松弛问题，设x*为可行解。
（2） 若x*是整数解，则算法停止。
（3） 若不是，增加一个线性不等式的限制条件使得所有的整数解满足，而x*不满足，返回第一步，重复直至停止。
设x*是松弛问题可行的基本解，且至少有一个基本变量。设N为非基本变量的集合。任意整数线性规划的解满足对任意x∈N，xi=0。这也是线性规划问题的解，一定与可行基本解x*相同。所以可行的整数解应该满足



∑j∈Nxj≥1




这是我们加入松弛问题的不等式。任何整数解都满足它，而x*不满足。
下面介绍具体的Gomory割平面法。
割平面法于1958年由高莫瑞（R.E.Gomory）首先提出，故又称Gomory割平面法。在割平面法中，每次增加的用于“切割”的线性约束称为割平面约束或Gomory约束。构造割平面约束的方法很多，但下面的方法是最常用的一种，它可以从相应线性规划的最终单纯形表中直接产生。
在松弛问题的最优单纯形表中，记Q为m个基变量的下标集合，K为n-m个非基变量的下标集合，则m个约束方程可表示为


xi+∑j∈Ka-ijxj=b-i,i∈Q（53）




而对应的最优解X*=(x*1，x*2,…,x*n)T，其中


x*j=b-j,j∈Q


0,j∈K（54）




若各b-j(j∈Q)皆为整数，则X*满足xj取整数(j=1,2,…,n)，因而就是纯整数规划的最优解； 若各b-j(j∈Q)不全为整数，则X*不满足xj取整数(j=1,2,…,n)，因而就不是纯整数规划的可行解，自然也不是原整数规划的最优解。
用割平面法（cutting plane approach）解整数规划时，若其松弛问题的最优解X*不满足xj取整数(j=1,2,…,n)，则从X*的非整分量中选取一个，用以构造一个线性约束条件，将其加入原松弛问题中，形成一个新的线性规划，然后求解。若新的最优解满足整数要求，则它就是整数规划的最优解； 否则，重复上述步骤，直到获得整数最优解为止。
为最终获得整数最优解，每次增加的线性约束条件应具备两个基本性质： ①已获得的不符合整数要求的线性规划最优解不满足该线性约束条件，从而不可能在以后的解中再出现； ②凡整数可行解均满足该线性约束条件，因而整数最优解始终被保留在每次形成的线性规划可行域中。
为此，若b-i0(i0∈Q)不是整数，在式（53）中对应的约束方程为



xi0+∑j∈Ka-i0,jxj=b-i0（55）




其中，xi0和xj(j∈K)按xj取整数(j=1,2,…,n)应为整数； b-i0按假设不是整数； a-i0,j(j∈K)可能是整数，也可能不是整数。
分解a-i0,j和b-i0成两部分： 一部分是不超过该数的最大整数，另一部分是余下的小数。即



a-i0,j=Ni0,j+fi0,j，Ni0,j≤a-i0,j且为整数，0≤fi0,j<1(j∈K)（56）

b-i0=Ni0+fi0，Ni0<b-i0且为整数，0<fi0<1（57）




把式（56）和式（57）代入式（55），移项得



xi0+∑j∈KNi0,jxj-Ni0=fi0-∑j∈Kfi0,jxj（58）




式（58）中，左边是一个整数，右边是一个小于1的数，因此有

fi0-∑j∈Kfi0,jxj≤0

即



∑j∈K(-fi0,j)xj≤-fi0（59）




现在考查线性约束条件（59）的性质。
一方面，由于式（59）中j∈K，所以，如将X*代入，各xj作为非基变量皆为0，将有



0≤-fi0




这和式（57）矛盾。由此可见，X*不满足式（59）。
另一方面，满足∑nj=1aijxj=bi(i=1,2,…,m)、xj≥0(j=1,2,…,n)和xj取整数(j=1,2,…,n)的任何一个整数可行解X一定也满足式（53）。式（55）是式（53）中的一个表达式，当然也满足。因而X必定满足式（58）和式（59）。由此可知，任何整数可行解一定能满足式（59）。
综上所述，线性约束条件（59）具备上述两个基本性质。将式（59）和max z=∑nj=1cjxj、∑nj=1aijxj=bi(i=1,2,…,m)、xj≥0(j=1,2,…,n)合并，构成一个新的线性规划。记R为原松弛问题可行域，R′为新的线性规划可行域。从几何意义上看，式（59）实际上对R做了一次“切割”，在留下的R′中，保留了整数规划的所有整数可行解，但不符合整数要求的X*被“切割”掉了。随着“切割”过程不断继续，整数规划最优解最终有机会成为某个线性规划可行域的顶点，作为该线性规划的最优解而被解得。
经验表明，实际解题时若从最优单纯形表中选择具有最小（分）数部分的非整分量所在行构造割平面约束，往往可以提高“切割”效果，减少“切割”次数。
5.6Gomory割平面法的MATLAB实现
Gomory割平面法的代码如下： 



%Gomory割平面法

function ［intx,intf］=Gomory(A,c,b,base)



%约束矩阵:A






%目标函数系数向量:c

%约束右端向量:b

%初始基向量:base

%目标函数取最小化时的自变量值:x

%目标函数的最小值:minf



sz=size(A);

nVia=sz(2);

n=sz(1);

xx=1:nVia;



if length(base)~=n

disp('基变量的个数要与约束矩阵的行数相等!');

mx=NaN;

mf=NaN;

return;

end



M=0;

sigma=-［transpose(c) zeros(1,(nVia-length(c)))］;

xb=b;

while 1

［maxs,ind］=max(sigma);

if maxs<=0

vr=find(c~=0,1,'last');

for l=1:vr

ele=find(base==l,1);

if(isempty(ele))

mx(l)=0;

else

mx(l)=xb(ele);

end

end

if max(abs(round(mx)-mx))<1.0e-7

intx=mx;

intf=mx*c;

return;

else

sz=size(A);

sr=sz(1);

sc=sz(2);

［max_x,index_x］=max(abs(round(mx)-mx));

［isB,num］=find(index_x==base);

fi=xb(num)-floor(xb(num));

for i=1:(index_x-1)






Atmp(1,i)=A(num,i)-floor(A(num,i));

end

for i=(index_x+1):sc

Atmp(1,i)=A(num,i)-floor(A(num,i));

end



%构建对单纯形法的初始表格

Atmp(1,index_x)=0;

A=［A zeros(sr,1);-Atmp(1,:) 1］;

xb=［xb;-fi］;

base=［base sc+1］;

sigma=［sigma 0］;



%对偶单纯形法迭代过程

while 1

if(xb)>=0

if max(abs(round(xb)-xb))<1.0e-7

%用对偶单纯形法求得了整数解

vr=find(c~=0 ,1,'last');

for l=1:vr

ele =find (base==l,1);

if(isempty(ele))

mx_1(l)=0;

else

mx_1(l)=xb(ele);

end

end

intx=mx_1;

intf=mx_1*c;

return;

else

sz=size(A);

sr=sz(1);

sc=sz(2);

［max_x,index_x］=max(abs(round(mx_1)-mx_1));

［isB,num］=find(index_x==base);

fi=xb(num)-floor(xb (num));

for i=1:(index_x-1)


Atmp(1,i)=a(num,i)-floor(A(num,i));



end

for i=(index_x+1):sc

Atmp(1,i)=a(num,i)-floor(a(num,i));

end



%下一次对偶单纯形法迭代的初始表格






Atmp(1,index_x)=0;

A=［A zeros(sr,1);-Atmp(1,:) 1］;

xb=［xb;-fi］;

base=［base sc+1］;

sigma=［sigma 0］;

continue;

end





%对偶单纯形法的换基变量过程

else

minb_1=inf;

chagB_1=inf;

sA=size(A);

［br,idb］=min(xb);

for j=1:sA(2)

if A(idb,j)<0

bm=sigma(j)/A(idb,j);

if bm<minb_1

minb_1=bm;

chagB_1=j;

end

end

end

sigma=sigma-A(idb,:)*minb_1;

xb(idb)=xb(idb)/A(idb,chagB_1);

A(idb,:)=A(idb,:)/A(idb,chagB_1);

for i=1:sA(1)

if i ~=idb

xb(i)=xb(i)-A(i,chagB_1)*xb(idb);

A(i,:)=A(i,:)-A(i,chagB_1)*A(idb,:);

end

end

base=chagB_1;

end

end

end

else

minb=inf;

chagB=inf;

for j=1:n

if A(j,ind)>0

bz=xb(j)/A(j,ind);

if bz <minb

minb=bz;






chagB=j;

end

end

end

sigma=sigma-A(chagB,:)*maxs/A(chagB,ind);

xb(chagB)=xb(chagB)/A(chagB,ind);

A(chagB,:)=A(chagB,:)/A(chagB,ind);

for i=1:n

if i~=chagB

xb(i)=xb(i)-A(i,ind)*xb(chagB);

A(i,:)=A(i,:)-A(i,ind)*A(chagB,:);

end

end

base(chagB)=ind;

end



M=M+1;

if (M==1000000)

disp('找不到最优解!');

mx=NaN;

minf=NaN;

return;

end

end





下面给出MATLAB解法的例子。
例58


min f(x)=x1-x2




使得



-x1+2x2≤2
2x1+x2≤4
x1,x2≥0,且x1,x2为整数




解： 首先引入人工变量x3、x4将约束条件转换为等式形式,即



min f(x)=x1-x2




使得



-x1+2x2+x3=2
2x1+x2+x4=4
x1,x2,x3,x4≥0,且x1,x2为整数




在MATLAB命令行输入下列命令： 



A=［-1210 ;2101］;

c=［1;-1］;






b=［2;4］;

［intx,intf］=Gomory(A,c,b,［3 4］)





结果如下： 



intx =

01

intf =

-1





例59用割平面法求解纯整数规划



max z=3x1-x2




使得



3x1-2x2≤3
5x1+4x2≥10
2x1+x2≤5
x1,x2≥0
x1,x2为整数




解： 引入松弛变量x3、x4、x5，将问题转换为标准形式，用单纯形法解其松弛问题，得最优单纯形表，见表54。


表54最优单纯形表



cj3－1000
CBXBbx1x2x3x4x5
3x113/7101/702/7
－1x29/701－2/703/7
0x431/700－3/7122/7
cj-zj00－5/70－3/7

由于b列各分数中x1=13/7有最大小数部分6/7，故从表54中第一行产生割平面约束。按照式（59），割平面约束为



-17x3-27x5≤-67（510a）




引入松弛变量x6，得割平面方程



-17x3-27x5+x6=-67（510b）




将式（510b）并入表54，然后用对偶单纯形法求解，得表55。


表55表54中第一行产生割平面约束后的新表




cj3－10000
CBXBbx1x2x3x4x5x6
3x113/7101/702/70
－1x29/701－2/703/70
0x431/700－3/7122/70
0x6－6/700－1/70［－2/7］1
cj-zj00－5/70－3/70
︙
3x11100001
－1x25/4010－1/40－5/4
0x35/2001－1/20－11/2
0x57/40001/41－3/4
cj-zj000－1/40－17/4


类似地，从表55中最后一个单纯形表的第四行产生割平面约束



-14x4-14x6≤-34（510c）




引入松弛变量x7，得割平面方程



-14x4-14x6+x7=-34（510d）




将式（510d）并入表55中最后一个单纯形表，然后用对偶单纯形法解之，得表56。


表56表55中第四行产生割平面约束后的新表




cj3－100000
CBXBbx1x2x3x4x5x6x7
3x111000010
－1x2201000－1－1
0x3400100－5－2
0x5100001－11
0x43000101－4
cj-zj－100000－4－1

表56给出的最优解(x1,x2,x3,x4,x5,x6,x7)T=(1,2,4,3,1,0,0)T已满足整数要求，故原整数规划问题的最优解为



x1=1，x2=2，max z=1




如果在先后构造的割平面约束式（510a）和式（510c）中，将各变量用原整数规划的决策变量x1和x2表示，则式（510a）和式（510c）


图54“切割”的几何意义示意图

成为x1≤1和x1+x2≥3。在这种形式下，“切割”的几何意义是显而易见的，如图54所示。
用MATLAB求解。先将限制条件转换为等式



max z=3x1-x2




使得



3x1-2x2+x3=3
5x1+4x2+x4=10
2x1+x2+x5=5
x1,x2,x3,x4,x5≥0




再如例58中，运用Gomory函数即可得到结果。
5.7分支定界法
分支定界法（branch and bound method）是一种隐枚举法（implicit enumeration）或部分枚举法，它不是一种有效算法，是在枚举法基础上的改进。分支定界法的关键是分支和定界。
若整数规划的松弛问题的最优解不符合整数要求，假设xi=b-i不符合整数要求，b-i是不超过b-i的最大整数，则构造两个约束条件： xi≤［b-i］和xi≥［b-i］+1。分别将其并入上述松弛问题中，形成两个分支，即两个后继问题。两个后继问题的可行域中包含原整数规划问题的所有可行解。而在原松弛问题可行域中，满足［b-i］<xi<［b-i］+1的一部分区域在以后的求解过程中被遗弃了，但它不包含整数规划的任何可行解。根据需要，各后继问题可以类似地产生自己的分支，即自己的后继问题。如此不断继续，直到获得整数规划的最优解。这就是所谓的“分支”。
所谓“定界”，是在分支过程中，若某个后继问题恰巧获得整数规划问题的一个可行解，则其目标函数值就是一个“界限”，可作为衡量处理其他分支的一个依据。因为整数规划问题的可行解集是它的松弛问题可行解集的一个子集，前者最优解的目标函数值不会优于后者最优解的目标函数值。所以，对于那些相应松弛问题最优解的目标函数值劣于上述“界限”值的后继问题，就可以剔除不再考虑了。当然，如果在以后的分支过程中出现了更好的“界限”，则以它来取代原来的界限，这样可以提高求解的效率。
“分支”为整数规划最优解的出现缩减了搜索范围，而“定界”则可以提高搜索的效率。经验表明，在可能的情况下，根据对实际问题的了解，事先选择一个合理的“界限”，可以提高分支定界法的搜索效率。
下面通过例子来阐明分支定界法的基本思想和一般步骤。
例510求解



max z=x1+x2




使得



x1+914x2≤5114
-2x1+x2≤13
x1,x2≥0
x1,x2取整数




解： 记整数规划问题为（IP），它的松弛问题为（LP）。图55中S为（LP）的可行域，黑点表示（IP）的可行解。用单纯形法解（LP），最优解为x1=3/2，x2=10/3，即点A，maxz=29/6。
（LP）的最优解不符合整数要求，可任选一个变量，如选择x1=3/2进行分支。由于最接近3/2的整数是1和2，因而可以构造两个约束条件



x1≥2（511a）




和



x1≤1（511b）




将式（511a）和式（511b）分别并入例510的松弛问题（LP）中，形成两个分支，即后继问题（LP1）和（LP2），分别由（LP）及式（511a）和（LP）及式（511b）组成。图56中S1和S2分别为（LP1）和（LP2）的可行域。不连通的域S1∪S2中包含了（IP）的所有可行解，S中被舍去的一部分S＼S1∪S2中不包含（IP）的任何可行解。



图55（LP）的可行域和（IP）的可行解




图56（LP1）和（LP2）的可行域



解（LP1），最优解为x1=2，x2=23/9，即点B，max z=41/9。点B仍不符合整数要求，再解（LP2）。（LP2）最优解为x1=1，x2=7/3，即点C，max z=10/3。点C也不符合整数要求，必须继续分支。
由于41/9>10/3，所以优先选择S1分支。因B点x1=2，而x2=23/9不符合整数要求，故可以构造两个约束条件



x2≥3（511c）




和



x2≤2（511d）




将式（511c）和式（511d）分别并入（LP1），形成两个新分支，即（LP1）的后继问题（LP11）和（LP12），分别由（LP1）及式（511c）和（LP1）及式（511d）组成。图57中S12为（LP12）的可行域。由于式（511c）和（LP1）不相容，故（LP11）无可行解，也就是说，（LP11）的可行域S11为空集，所以只需考虑后继问题（LP12）。
在S12上解（LP12），最优解为x1=33/14，x2=2，即图57中点D，max z=61/14。
对于原整数规划（IP）来说，至此还剩两个分支： 后继问题（LP2）和（LP12）。因为（LP12）的最优解目标函数值比（LP2）大，所以优先考虑对（LP12）进行分支。
两个新约束条件为



x1≥3（511e）




和



x1≤2（511f）




类似地，形成（LP12）的两个后继问题（LP121）和（LP122）。图58中S121和S122分别为它们的可行域，其中S122是一条直线段。
（LP121）的最优解是x1=3，x2=1，即图58中的点E，max z=4； （LP122）的最优解是x1=2，x2=2，即图58中的点F，max z=4。这两个解都是（IP）的可行解，且目标函数值相等。至此，可以肯定两点： 第一，在S121和S122中不可能存在比E点和F点更好的（IP）的可行解，因此不必再在它们中继续搜索； 第二，既然点E和点F都是（IP）的可行解，它们的目标函数值z=4就可看作（IP）最优解的目标函数值的一个界限（对于最大化问题，是下界； 对于最小化问题，是上界）。



图57（LP12）的可行域




图58（LP121）和（LP122）的可行域



现在，尚未检查的后继问题只有（LP2）了。但（LP2）的最优解的目标函数值是10/3，比界限4小。因此，S2中不存在目标函数值比4大的（IP）的可行解，也就是说，不必再对（LP2）进行分支搜索了。
综上所述，我们已经求得了整数规划（IP）的两个最优解。它们分别是x1=3，x2=1和x1=2，x2=2，max z=4。
上述分支定界法求解的过程可用图59表示。


图59分支定界法求解的过程



分支定界法解整数规划的一般步骤如下。
步骤1称整数规划问题为问题A，它的松弛问题为问题B，以zb表示问题A的目标函数的初始界（如已知问题A的一个可行解，则可取它的目标函数值为zb）。对最大化问题A，zb为下界； 对最小化问题A，zb为上界。解问题B。转步骤2。
步骤2如问题B无可行解，则问题A也无可行解； 如问题B的最优解符合问题A的整数要求，则它就是问题A的最优解。对于这两种情况，求解过程到此结束。如问题B的最优解存在，但不符合问题A的整数要求，则转步骤3。
步骤3对问题B，任选一个不符合整数要求的变量进行分支。设选择xj=b-j，且设b-j为不超过b-j的最大整数。对问题B分别增加下面两个约束条件中的一个： 



xj≤［b-j］和xj≥［b-j］+1


从而形成两个后继问题。解这两个后继问题。转步骤4。
步骤4考查所有后继问题，如其中有某几个存在最优解，且其最优解满足问题A的整数要求，则以它们中最优的目标函数值和界zb作比较。若比界zb更优，则以其取代原来的界zb，并称相应的后继问题为问题C。否则，原来的界zb不变。转步骤5。
步骤5不属于C的后继问题中，称存在最优解且其目标函数值比界zb更优的后继问题为待检查的后继问题。
若不存在待检查的后继问题，当问题C存在时，问题C的最优解就是问题A的最优解； 当问题C不存在时，和界zb对应的可行解就是问题A的最优解。zb即为问题A的最优解的目标函数值，求解到此结束。
若存在待检查的后继问题，则选择其中目标函数值最优的一个后继问题，改称其为问题B。回到步骤3。
分支定界法是求解整数规划的较好方法，很多求解整数规划的计算机软件是根据分支定界法原理编写的，同时这种方法也适用于求解混合整数规划问题，在实际中有着广泛应用。
5.8分支定界法的MATLAB实现
根据分支定界算法的步骤，可以写出其具体实现代码： 



% By Sherif A. Tawfik, Faculty of Engineering, Cairo University

% ［x,val,status］=IP1(f,A,b,Aeq,beq,lb,ub,M,e)



% this function solves the following mixed-integer linear programming problem

%min f*x

%subject to

%A*x <=b

%Aeq * x = beq

%lb <= x <= ub

%M is a vector of indeces for the variables that are constrained to be integers

%e is the integarilty tolerance



% the return variables are :

% x : the solution

% val: value of the objective function at the optimal solution

% status =1 if successful

%=0 if maximum number of iterations reached in he linprog function

%=-1 if there is no solution



% Example:

%maximize 17 x1 + 12 x2

%subject to

%10 x1 + 7 x2 <=40

%x1 +x2 <= 5

%x1, x2 >=0 and are integers

% f=［-17, -12］; %take the negative for maximization problems

% A=［ 107; 1 1］;

% B=［40; 5］;

% lb=［0 0］;

% ub=［inf inf］;

% M=［1,2］;

% e=2^-24;






% ［x v s］= IP(f,A,B,［］,［］,lb,ub,M,e)





function ［x,val,status］=intprog(f,A,b,Aeq,beq,lb,ub,M,e)

options = optimset('display','off');

bound=inf; % the initial bound is set to +ve infinity

［x0,val0］=linprog(f,A,b,Aeq,beq,lb,ub,［］,options);

［x,val,status,b］=rec(f,A,b,Aeq,beq,lb,ub,x0,val0,M,e,bound); % a recursive function that processes the BB tree



function ［xx,val,status,bb］=rec(f,A,b,Aeq,beq,lb,ub,x,v,M,e,bound)

options = optimset('display','off');

% x is an initial solution and v is the corressponding objective function value

% solve the corresponding LP model with the integarily constraints removed



［x0,val0,status0］=linprog(f,A,b,Aeq,beq,lb,ub,［］,options);



% if the solution is not feasible or the value of the objective function is

% higher than the current bound return with the input intial solution



if status0<=0 | val0 > bound

xx=x; val=v; status=status0; bb=bound;

return;

end





% if the integer-constraint variables turned to be integers within the

% input tolerance return



ind=find( abs(x0(M)-round(x0(M)))>e );

if isempty(ind)

status=1;

if val0 < bound% this solution is better than the current solution hence replace

x0(M)=round(x0(M));

xx=x0;

val=val0;

bb=val0;

else

xx=x;% return the input solution

val=v;

bb=bound;

end

return

end



% if we come here this means that the solution of the LP relaxation is






% feasible and gives a less value than the current bound but some of the

% integer-constraint variables are not integers.

% Therefore we pick the first one that is not integer and form two LP problems

% and solve them recursively by calling the same function (branching)





% first LP problem with the added constraint that Xi <= floor(Xi) , i=ind(1)



br_var=M(ind(1));

br_value=x(br_var);

if isempty(A)

［r c］=size(Aeq);

else

［r c］=size(A);

end

A1=［A ; zeros(1,c)］;

A1(end,br_var)=1;

b1=［b;floor(br_value)］;



% second LP problem with the added constraint that Xi >= ceil(Xi) , i=ind(1)



A2=［A ;zeros(1,c)］;

A2(end,br_var)=-1;

b2=［b; -ceil(br_value)］;





% solve the first LP problem



［x1,val1,status1,bound1］=rec(f,A1,b1,Aeq,beq,lb,ub,x0,val0,M,e,bound);

status=status1;

if status1 >0 & bound1<bound % if the solution was successfull and gives

a better bound

xx=x1;

val=val1;

bound=bound1;

bb=bound1;

else

xx=x0;

val=val0;

bb=bound;

end



% solve the second LP problem



［x2,val2,status2,bound2］=rec(f,A2,b2,Aeq,beq,lb,ub,x0,val0,M,e,bound);



if status2 >0 & bound2<bound % if the solution was successfull and gives






a better bound.

status=status2;

xx=x2;

val=val2;

bb=bound2;

end




例511


min z=-5x1-x2




使得



3x1+x2≥9
x1+x2≥5
x1+8x2≥8
x1,x2≥0且为整数




解： 在MATLAB输入命令



f=［-5,-1］;

A=［3,1;1,1;1,8］;

b=［9;5;8］;

lb=［0 0］;

ub=［inf, inf］;

M=［1,2］;

e=2^-24;

［x,v,s］=intprog(f,A,b,［］,［］,lb,ub,M,e)





得到



x =

3

0

v =

-15.0000

s =

1





故最小值为15，此时x1=3，x2=0。
例512


max z=x1+x2




使得



x1+914x2≤5114
-2x1+x2≤13
x1,x2≥0
x1,x2取整数




解： 当求最大值时，先添加负号，相当于求



min z=-x1-x2




这样即可化为求解的标准型。
在MATLAB中输入下列命令： 



f=［-1,-1］;

A=［1,9/14;-2,1］;

b=［51/14;1/3］;

lb=［0 0］;

ub=［inf, inf］;

M=［1,2］;

e=2^-24;

［x,v,s］=intprog(f,A,b,［］,［］,lb,ub,M,e)





结果为



x =

3

1

v =

-4.0000

s =

1





故最大值为4。
5.9整数规划的解法
01型整数规划是一种特殊的整数规划，若含有n个变量，则可以产生2n个可能的变量组合。当n较大时，采用完全枚举法解题几乎是不可能的。已有的求解01型整数规划的方法一般都属于隐枚举法。
在2n个可能的变量组合中，往往只有一部分是可行解。只要发现某个变量组合不满足其中一个约束条件时，就不必再去检验其他约束条件是否可行。对于可行解，其目标函数值也有优劣之分。若已发现一个可行解，则根据它的目标函数值可以产生一个过滤条件（filtering constraint），即对于目标函数值比它差的变量组合就不必再去检验它的可行性。在以后的求解过程中，每次发现比原来更好的可行解，即以此替换原来的过滤条件。上述做法都可以减少运算次数，使最优解能较快地被发现。
例513求解01整数规划



max z=3x1-2x2+5x3




使得


x1+2x2-x3≤2
x1+4x2+x3≤4
x1+x2≤3
4x2+x3≤6
x1,x2,x3=0或1（512）




解： 求解过程可以用表表示（见表57）。


表57例513的求解过程




x1,x2,x3z值约束条件

abcd过滤条件
(0,0,0)0√√√√z≥0
(0,0,1)5√√√√z≥5
(0,1,0)－2
(0,1,1)3
(1,0,0)3
(1,0,1)8√√√√z≥8
(1,1,0)1
(1,1,1)6

所以，最优解(x1,x2,x3)T=(1,0,1)T,max z=8。
采用上述算法，实际只做了20次运算。
为了进一步减少运算量，常按目标函数中各变量系数的大小顺序重新排列各变量，以使最优解有可能较早出现。对于最大化问题，可按由小到大的顺序排列； 对于最小化问题，则相反。为此例513可写成下列形式:



max z=5x3+3x1-2x2




使得


-x3+x1+2x2≤2
x3+x1+4x2≤4
x1+x2≤3
x3+4x2≤6
x3,x1,x2=0或1（513）




求解时先令排在前面的变量取值为1，如本例中可取(x3,x1,x2)=(1,0,0)，若不满足约束条件时，可调整取值为（0,1,0）； 若仍不满足约束条件，可退为取值（0,0,1）等，依此类推。据此改写后模型的求解过程可见表58。


表58改写后模型的求解过程



(x3,x1,x2)z值约束条件

abcd过滤条件
(0,0,0)0√√√√z≥0
(1,0,0)5√√√√z≥5
(1,1,0)8√√√√z≥8

从目标函数方程看到，z值已不可能再增大，(x3,x1,x2)=(1,1,0)即为本例的最优解。
采取这样的形式用上述方法解此例，可以很大程度减少运算次数。一般问题的规模越大，这样做的好处就越明显。
5.1001整数规划的MATLAB实现
对于混合整数规划问题，MATLAB Optimization Toolbox中给出了INTLINPROG函数来解决。可以看到，INTLINPROG可以解决以下（混合）整数线性规划问题。



function ［x, fval, exitflag, output］ = intlinprog(f, intcon, A, b, Aeq, beq, lb, ub, x0, options)

%INTLINPROG Mixed integer linear programming.

%

%X = INTLINPROG(f,intcon,A,b) attempts to solve problems of the form

%

%min f'*xsubject to:A*x<= b

%xAeq*x = beq

%lb <= x <= ub

%x(i) integer, where i is in the index

%vector intcon (integer constraints)




函数INTLINPROG能够求



min fTx




满足条件



x(在intcon中的位置上的值)是整数

Ax≤b

Aeq*x=beq

lb≤x≤ub




在要求01整数规划问题时，可令intcon包含所有x的元素，并且lb和ub分别是全为0、1的向量即可。当然，该函数对于所有的（混合）整数线性规划都是适用的，应用十分广泛。
更多该函数的详细代码可以在intlinprog.m中查阅，不在此赘述。
下面是一个使用INTLINPROG解01整数规划的例子。
例514


min f(x)=x1+2x2+3x3+x4+x5




使得



2x1+3x2+5x3+4x4+7x5≥8
x1+2x2+4x3+2x4+2x5≥5
x1,x2,x3,x4,x5=0或1





在MATLAB中输入下列命令： 



f=［1;2;3;1;1］;

intcon=1:5;

A=［-2,-3,-5,-4,-7;-1,-2,-4,-2,-2］;

b=［-8;-5］;

lb=［0,0,0,0,0］;

ub=［1,1,1,1,1］;

［x,fval］ = intlinprog(f,intcon,A,b,［］,［］,lb,ub)





得到



x =

1

0

0

1

1

fval =



3





5.11整数规划解决旅行商问题的MATLAB实例

这个例子将展示如何使用二进制整数规划来解决经典的旅行商问题。这个问题涉及通过一组城市（200个城市）找到最短的封闭旅游路径。
其中包含要解决的问题如下： 对于所有不同的城市站点，生成所有可能的行程； 计算每次旅行的距离； 最小化成本函数即每次旅行距离之和。
将问题公式化，设置二进制的决策变量，并与每个行程相关联： 1表示在行程中的路径，0表示不在行程中的路径。为了确保行程经过每个站点，设置线性限制条件使得每个站点链接两条路径，即一条离开一条到达的路径。
（1） 生成停止条件。
将这组城市的连线粗略地看作多边形，在其中设置随机停止： 



load('usborder.mat','x','y','xx','yy');

rng(3,'twister') % Makes a plot with stops in Maine & Florida, and is reproducible

nStops = 200; % You can use any number, but the problem size scales as N^2

stopsLon = zeros(nStops,1); % Allocate x-coordinates of nStops

stopsLat = stopsLon; % Allocate y-coordinates

n = 1;

while (n <= nStops)

xp = rand*1.5;

yp = rand;

if inpolygon(xp,yp,x,y) % Test if inside the border

stopsLon(n) = xp;

stopsLat(n) = yp;

n = n+1;

end

end





（2） 计算两个站点之间的距离。
生成所有的路径： 



idxs = nchoosek(1:nStops,2);





因为有200个站点，所以一共有19900条路径，即有19900个二进制变量(# variables = 200 choose 2)。
计算所有的旅行距离，假设地球是平的，以便使用毕达哥拉斯法则。



dist = hypot(stopsLat(idxs(:,1)) - stopsLat(idxs(:,2)), ...

stopsLon(idxs(:,1)) - stopsLon(idxs(:,2)));

lendist = length(dist);




（3） 创建图形并绘制地图。
绘图表示问题（如图510所示），其中点为站点，边为路径。



G = graph(idxs(:,1),idxs(:,2));

figure

hGraph = plot(G,'XData',stopsLon,'YData',stopsLat,'LineStyle','none','NodeLabel',{});






hold on

% Draw the outside border

plot(x,y,'r-')

hold off







图510旅行商问题城市地图






①
1英里=1.609千米。



（4） 创建约束条件。
创建线性约束条件，即每个站点都有两个关联的路径，因为每个站点必须有一个到达路径和一个离开站点的路径。



Aeq = spalloc(nStops,length(idxs),nStops*(nStops-1)); % Allocate a sparse matrix

for ii = 1:nStops

whichIdxs = (idxs == ii); % Find the trips that include stop ii

whichIdxs = sparse(sum(whichIdxs,2)); % Include trips where ii is at either end

Aeq(ii,:) = whichIdxs'; % Include in the constraint matrix

end

beq = 2*ones(nStops,1);





增加二进制限制，将intcon参数设置为决策变量的个数，每个变量下界为0，上界为1。



intcon = 1:lendist;

lb = zeros(lendist,1);

ub = ones(lendist,1);





（5） 用intlinprog进行优化，将结果可视化。
创建一个用来求解路径的新图。为此，在某些值不是整数的情况下对解决方案进行舍入，并将结果值转换为逻辑值。



opts = optimoptions('intlinprog','Display','off');

［x_tsp,costopt,exitflag,output］ = intlinprog(dist,intcon,［］,［］,Aeq,beq,lb,ub,opts);

x_tsp = logical(round(x_tsp));

Gsol = graph(idxs(x_tsp,1),idxs(x_tsp,2));

hold on

highlight(hGraph,Gsol,'LineStyle','-')

title('Solution with Subtours')





如图511所示，该解决方案有几个子环路。目前为止指定的约束并不能阻止这些子环路的发生。为了防止任何可能的子环路发生，需要大量的不等式约束。


图511旅行商最佳旅游线路



（6） 避免子环路的限制条件。
由于无法添加所有的子环路的约束，这里采用迭代方法。检测当前解决方案中的子环路，然后添加不等式约束以防止这些特定的子环路发生。这样就可以在几个迭代中找到一个合适的路径。
消除带有不等式约束的子函数。一个这样做的例子是，如果在一个子环路中有5个点，那么有5条线连接这些点来创建子环路。通过实现一个不等式约束来消除子环路，即这5个点之间必须有小于或等于4条线。
更重要的是，找到这5个点之间的所有直线，并增加约束条件，使这些直线不超过4条。如果一个解决方案中存在5条或更多的行，那么该解决方案将具有一个子环路（具有n个节点和n条边的图总是包含一个循环），如图512所示。
加入线性不等式约束以消除子环路，并重复调用，直到剩下一个环路。



tourIdxs = conncomp(Gsol);

numtours = max(tourIdxs); % number of subtours

fprintf('# of subtours: %d＼n',numtours);



A = spalloc(0,lendist,0); % Allocate a sparse linear inequality constraint matrix

b = ［］;

while numtours > 1 % Repeat until there is just one subtour

% Add the subtour constraints

b = ［b;zeros(numtours,1)］; % allocate b

A = ［A;spalloc(numtours,lendist,nStops)］; % A guess at how many nonzeros to allocate

for ii = 1:numtours

rowIdx = size(A,1) + 1; % Counter for indexing

subTourIdx = find(tourIdxs == ii); % Extract the current subtour

%The next lines find all of the variables associated with the

%particular subtour, then add an inequality constraint to prohibit

%that subtour and all subtours that use those stops.

variations = nchoosek(1:length(subTourIdx),2);

for jj = 1:length(variations)

whichVar = (sum(idxs==subTourIdx(variations(jj,1)),2)) & ...

(sum(idxs==subTourIdx(variations(jj,2)),2));

A(rowIdx,whichVar) = 1;

end

b(rowIdx) = length(subTourIdx) - 1; % One less trip than subtour stops

end



% Try to optimize again

［x_tsp,costopt,exitflag,output］ = intlinprog(dist,intcon,A,b,Aeq,beq,lb,ub,opts);

x_tsp = logical(round(x_tsp));

Gsol = graph(idxs(x_tsp,1),idxs(x_tsp,2));



% Visualize result

hGraph.LineStyle = 'none'; % Remove the previous highlighted path

highlight(hGraph,Gsol,'LineStyle','-')

drawnow



% How many subtours this time?

tourIdxs = conncomp(Gsol);

numtours = max(tourIdxs); % number of subtours

fprintf('# of subtours: %d＼n',numtours)

end



title('Solution with Subtours Eliminated');

hold off






图512旅行商避免子环路后的最佳旅游线路