深入浅出列生成算法
本文尽量避免数学公式,使用文字解释列生成算法的原理,争取让读者能形成直观上的理解。
为什么需要了解列生成算法的原理
- 列生成算法无法简单地调用第三方库来使用,必须根据具体问题,构造不同的算法模型。
- 只有了解了原理,才能在踩到各种坑时,有所针对地去优化各种细节。不然只能抓瞎或者抓腮。
列生成算法原理
列生成算法可以从两个视角来理解:对偶角度和单纯形算法角度。
对偶角度
啥是对偶
这里简单过一下对偶的概念。
假设有个长得很标准的线性规划问题:
那么,它的对偶问题为:
下面我们都以这个问题来讨论,即说到原问题时,默认是一个最小化问题;说到对偶问题时,默认是一个最大化问题。
怎么理解这个对偶关系呢?借用经济学方面的话来说,假设原问题的目标是让成本最小,那么对偶就是让收入最大。更确切地讲,是:
- 原问题丶:保证收入不低于某个值的条件下,使成本最小化。
- 对偶问题:保证成本不高于某个值的条件下,使收入最大化。
那个丶纯粹是为了对齐,忽略之……
可以看到,原问题和对偶问题其实就是一个问题:目标净收益最大。只是一个是约束收入优化成本,一个是约束成本优化收入。角度不同而已。体现在公式上,就是原问题的变量对应对偶问题的约束,目标系数对应约束边界,约束矩阵倒转过来。
另外,关于对偶,一个比较重要的特性是:原问题的最优值与对偶问题的最优值相等。
从对偶角度看列生成算法
列生成算法主要用途在于求解变量多,但是大部分变量将会取值为0的线性规划问题。总体思路是先忽略大部分变量,构造一个只使用小部分变量的模型(其余变量相当于值为0),这样就能很快求出一个解。然后寻找模型外的变量,找到能够让目标值更优的变量,加入模型再次求解。重复这个过程直至找不到更好的变量。
这个过程的关键问题在于,怎么评估模型外的变量是否能让目标值更优。
我们从对偶的角度来研究这个问题。
原问题的变量对应对偶问题的约束。所以原问题新增变量,相当于对偶问题新增约束。
原问题新增变量 -> 对偶问题新增约束
由于对偶问题是个最大化问题,所以对偶问题新增约束后,显然最优值不变或变差,也就是不变或变小。从常理上看,约束越多,最优值越差嘛。
而前面提到的,原问题的最优值等于对偶问题的最优值。也就是说,如果对偶问题最优值不变,那么原问题最优值也不变;如果对偶问题最优值变小,那么原问题最优值也变小。而我们需要的正是让原问题的最优值变小。
所以问题变为如何尽量避免新增的约束没有改变最优值。设想一下,当加入新约束时,如果当前对偶的最优解没有违反新的约束,那么这个解仍然会是新增约束后的对偶问题的最优解,最优值将不变。
因此,我们要找的新增的约束,要和当前最优解冲突。
整条逻辑链为:
新增变量后原问题最优解变小 -> 新增约束后对偶问题最优解变小 -> 新增约束前的最优解不在新增约束后的可行域 -> 新增约束前的最优解不满足新增的约束
一行对偶问题的约束的公式为:
假设最优解为w*,那么违反约束的条件为:
变换一下,变成:
左侧的式子,叫做的reduced cost,也叫做检验数。
通过分析,我们知道,只要加入reduced cost小于0的对偶约束(从而加入了原问题对应的变量)即可。
很自然的想法是,我们更倾向于找到reduced cost最小的一个或几个变量加入,也就是最好能找到最小化reduced cost的新约束:
这里就出现了一个新的最优化问题。这个问题叫做列生成的子问题(sub problem)。其中w*是已知的,未知量是c和a。c和a是和问题的应用场景有关的,需要根据实际场景来构造c和a的约束条件。所以子问题无法通用地求解,只能根据具体问题选择不同的方法求解。
当所有未加入模型的变量的reduced cost都大于等于0时,目标值无法再优化,说明我们已得到最优解。
另外,熟悉对偶问题经济学含义的同学会知道,reduced cost是指产品的差额成本。那么显然要新增的是差额成本为负的产品了。这是另一种理解列生成的思路。
单纯形算法角度
对偶角度给出了一个偏感性的方式来理解列生成算法。换个视角,从单纯形算法角度上看,则是单纯形算法本身,为了更高效地求解包含大量变量的问题,自然地扩展为列生成算法。
相信有不少人被单纯形算法虐得有心理阴影——公式复杂,手工计算量也巨大……
其实,如果我们先不看细节,单纯形的核心原理并没有那么难以理解。下面讲解时不会很严谨,理解算法框架就够了。严谨的过程请参阅运筹学相关书籍。
单纯形算法
众所周知,单纯形算法有一个几何上的解释:
- 线性规划是一个凸优化问题,局部最优解就是全局最优解。
- 线性规划的解空间是一个n维的凸多面体。最优解在这个凸多面体的某个顶点上。
- 单纯形算法从一个初始顶点开始,不断沿着邻边找更好的顶点。
- 当一个顶点四周没有更好的顶点时,这个顶点就是最优解。
整个过程就像水沿着一条蜿蜒的沟渠流下,最终汇聚到最低点一样。
问题是,这里面的几何概念和代数公式怎么对应?
这里用不严谨(但更容易理解)的语言说明一下:
- 边界:解空间是由不等式约束(包括变量非负这些约束)围起来的一块空间区域。当点p使得若干个不等式取等号时,那么点p就在约束边界的超平面上。这个边界可能是一个面、边、顶点。
- 顶点:顶点会让尽量多的约束取等号。也就是说,顶点是由若干个改为等号的约束组成的方程组的解。我们叫这个方程组为约束边界方程组。
- 沿着边:约束边界方程组去掉一个方程,其解集就变成与顶点邻接的一条边。再取一条原方程组外的约束条件加入,所得到的解就是相邻的顶点。简单说,就是约束边界方程组中替换掉一个方程,形成的新方程组解出来就是相邻的顶点。
这里涉及到通过让约束取等号来求边界的操作,而不等式乱糟糟地混在方程型的约束和变量非负约束里,会使这里的分析比较困难。所以使用单纯形算法之前,都会通过引入松弛变量、剩余变量和人工变量等方法(这一步在这里不重要,不详细展开了),将线性规划转换成如下标准形式:
标准形式中只有变量非负约束包含不等式,其他约束都是等式。这样我们就可以很容易地做边界相关的计算了。假设变量数量为n,等式约束数量为m。通过转换而来的标准形式都会有n > m。那么,我们知道,只要让n-m个变量等于0,剩下的m个变量就可以通过这m个等式联立方程组(约束边界方程组)求出一个解(简单起见,不考虑无解,无数解这些边缘条件)。这个解就是一个顶点。
这里约束边界方程组中的m个变量叫作基变量,固定值为0的n-m个变量叫作非基变量。
沿着边找相邻顶点,就是取一个被固定为0的非负约束,也就是一个非基变量(这个操作叫入基),替换掉一个基变量(这个操作叫出基,这个变量出基后就固定值为0),然后重新求解一个顶点。
入基操作需要选择入基变量,选择的依据是这个变量在目标函数中的下降速度,也就是这个变量增加1时,目标值减少多少。经过推导可知,下降速度的计算公式刚好是检验数(reduced cost)。这里就和对偶的视角联系起来了。
出基操作这里就不细说了,大致的思路是在约束条件下,旧的基变量有一部分会随着入基变量的增长而下降,其中最先下降到0的旧的基变量就会被选为出基变量。
整个单纯形算法的计算步骤是:
- 选取基变量和非基变量,简单能出初始解就好。
- 计算所有非基变量的reduced cost,找到最小且为负值的那个作为入基变量。如果reduced cost都大于等于0,迭代终止。
- 选出基变量
- 解约束边界方程组,回到步骤2
从单纯形算法角度看列生成算法
在单纯形的步骤2,需要计算所有非基变量的RC。找到最小的那个。当变量个数很多的时候,这一步就成为了算法运行时间瓶颈。
在一些情况下,通过巧妙构造问题,可以让这一步不需要遍历所有变量。甚至我们都不需要知道有多少变量,只要能在每次迭代的时候生成一个或者多个变量,提升优化效果就可以了。
由于不需要遍历所有变量,所以一开始就不需要使用所有变量,只需要使用一组能产生初始解的初始变量构成线性优化问题即可。这种只使用部分变量的模型被称为原问题的restricted master problem(RMP)。
每次迭代时,生成一个或多个让reduced cost最小的变量加入RMP。这个生成步骤就是求解子问题。不断加入新变量直到没有小于0的reduced cost的变量时就达到最优解。
到这里就和对偶角度分析的结果一致了。
下面是单纯形算法与列生成算法简要流程图的对比,可以看到,两者的结构是一样的。
一般来说,我们不会手搓单纯形算法,所以正常都是直接调用单纯形算法库解RMP,然后做列生成,再跑RMP,直到达到最优。
一个经典例子:Cutting Stock Problem
这是一个列生成算法的经典例子。
原纸卷每个长17m,顾客们分别需要25个3m长,20个5m长,18个7m长的纸卷。
问:如何切割使消耗的原纸卷数量最少?
令一个原纸卷的切割方案集合为:
P = {(a, b, c) | 3a + 5b + 7c <= 17}
其中,a是一个原纸卷切割出的3m纸卷数量,b是5m纸卷数量,c是7m纸卷数量。
我们用变量x(abc)表示使用切割方案(a, b, c)的原纸卷数量。
显然,一个变量与一个原纸卷切割方案一一对应。建模如下:
这里故意不适用传统的下标序号标记,意在突出我们不需要对变量编号,只需要知道变量在对应在什么集合上,如何通过集合中的元素生成变量就行了。
初始解很好找。比如说我们可以取25个原纸卷按照方案(1, 0, 0)切割,20个原纸卷按照方案(0, 1, 0)切割,18个原纸卷按照方案(0, 0, 1)切割。这当然会有很多浪费。但是初始解可行就可以了,浪费的部分会在下面的迭代中优化掉。
接下来要生成变量。变量与切割方案一一对应的。所以是要找出一个切割方案(a, b, c),使得reduced cost最小。
其中w1、w2和w3分别为约束R1、R2和R3的对偶值。
约束条件除了a、b、c非负外,还需要满足切割后的纸卷长度综合小于或者等于原纸卷的长度。
这样子问题就构造好了。求解子问题得到新增变量。然后迭代直到最优。具体计算这里不展开了。
整数规划求解
前面提到的单纯形算法和列生成算法求解的都是线性规划。在实际应用中,一般还会需要求解整数规划。也就是变量都约束为整数的线性规划。
这里先提一个概念:整数规划的线性松弛,整数规划问题,不考虑整数约束,剩下的约束条件和目标组成的线性规划问题。
其实我们并没有很好的方法直接求解整数规划,通常都是不断地调整并求解线性松弛,最后找到最优整数解。
分支定界
分支定界是一个用来求整数优化问题的框架。其实思路很简单,就是采用类似二分法的技巧,在线性解空间中暴力搜索整数解。
首先,求解线性松弛得到线性解。取一个变量x进行分支。比如x线性解值为1.2,那么产生 x <= 1 和 x >= 2 两个分支。将这两个条件分别加入到线性松弛,得到两个线性规划。再求解这两个线性规划,两个又分别分支……直到求得最优解。
有些情况下判断一个解是否为最优解是有方法的,所以不用搜索所有分支。但是,分支定界在最坏情况下的时间复杂度仍是指数级。为了防止运行时间过长,一般使用分支定界时还会额外加一些终止条件,比如回溯次数限制、运行时间限制、找到第一个整数解就结束等。
分支定价
分支定价就是在分支定界框架中,使用列生成算法来求解每个分支节点的算法。
不过这里,除了根节点,其他节点不用从头开始生成新变量,继承父节点用到的变量即可,这样可以节省很多重复生成变量的过程。
在01整数规划中,还有更简化的方法,每次列生成得到线性松弛的最优解后,找出值最接近1的变量,新加这些变量等于1的约束,继续跑列生成,直到找出所以值为1的变量。剩下的自然都是0了。这个方法可以加入回溯,也可以不回溯,出现无解就直接结束……据说不回溯也很少出现无解的情况。
一些相关问题
退化问题/类退化问题
通过RC找的新变量不一定能让目标值变得更好,仍然存在不变的可能。极小概率的情况下,单纯形算法可能会有入基变量和出基变量循环出现的情况。由于我们肯定是调用线性规划库来跑单纯形的,所以不用考虑这个……
列生成没有出基操作,不会出现循环。但是有一些改进会剔除冗余变量,这时就会有极小概率会出现循环了。这种情况不需要费心去处理,玄学调参降低出现概率,并设置最大迭代次数等强制终止条件,确保能终止就好。
最恶心的情况是没有循环,但是长时间没有提升目标值的情况。这其实是算法卡在一个拐点上了,只要过了拐点就能开始提升。特别是在一些约束较强的问题(比如密集的排班问题)中,使用某些启发式算法或者手工做出来的初始解就很容易出现这种情况。
而我们为了避免算法跑太久,通常会设置多次迭代没有提升就结束的条件。这可能使算法从拐点出发后,几次迭代无优化就直接结束了。
这种情况无法完美解决。简单的就是调参加更巧妙地设置结束条件,通过多次试验尽量让算法能跨过这个拐点。还有另一个技巧是可以适当地给约束边界加一下噪音,比如说小于等于1的约束,可以放宽到小于等于1.0001。这样从初始解出发迭代时,由于边界宽松了一些,变量可以有些许变化,会让目标值有一些微小的提升,帮助判断是否需要结束迭代。
CPLEX计算reduced cost的问题
使用CPLEX时,我们可以很方便的设置变量的上下界。比如设置 0 <= x <= 1。这时,x <= 1这部分是会影响reduced cost的值的。而CPLEX接口计算的是没有考虑这个条件的……所以可能你自己手搓代码出来reduced cost和CPLEX接口出来的reduced cost不相等。
更严重的是,可能你会忘了 x <= 1 这个条件,导致列生成的过程中算错reduced cost。
这个问题其实影响不大,主要是会干扰一些计算过程正确性的验算。
如何验证子问题有没有严重问题
没有做分支操作时,线性松弛的目标值如果变差,说明子问题可能出现了一些很蠢的问题。
线性规划求解算法选择
每次迭代求解线性规划时,选用不同的算法会影响求解时间。根据经验:
- 增加少量列时(列生成),使用单纯形算法(Primal)。
- 增加少量约束时(分支),使用对偶单纯形算法(Dual)。
- 其他情况,酌情使用对偶单纯形算法或者内点法。通过试验决定。
- 对偶单纯形算法快的时候很快,慢的时候很慢。
- 内点法速度比较稳定。
以上也并不是所有场景通用的。应当针对具体问题,反复试验来确定使用什么算法。
迭代中剔除冗余变量
Reduced cost可以用来评估变量的“有用”程度。越小表示变量越有用,越大表示变量越无关紧要。
列生成迭代次数较多后,变量数量会越来越多,从而每次迭代的运行速度越来越低。可以设定一个变量规模上限,当变量数量大于上限时,从模型中去掉reduced cost最大的那些变量。