[Python知识库]拉格朗日松弛求解整数规划浅析（附Python代码实例）

拉格朗日松弛是一种常见的求解大规模整数规划/混合整数规划的有效方法。目前网络上关于拉格朗日松弛的解读多数是从拉格朗日松弛求解非线性规划问题的角度着手，很少有关于拉格松弛求解整数规划的内容。另一方面在实践层面上拉格朗日松弛方法 有着诸多需要注意的问题，而教科书理论证明往往为了proof的漂亮很少讨论这些缺陷和问题。
我在博士阶段先后两次修了Nonlinear programming的课程，主讲老师是Peter Luh老师（陆宝森老师）。陆老师是研究拉格朗日松弛求解整数规划领域的前驱者，我十分有幸和陆老师讨论交流过很多问题，后边自己也动手实践过拉格朗日松弛算法，并发表了相关的paper，因此这里我把一些经验分享出来。

0 大规模整数规划/混合整数规划的 decomposition method 简介

由于大规模整数规划/混合整数规划往往是NP-hard的，所以我们很难直接对这些问题求解。每当遇到这些复杂问题的时候，一个基本思路就是对复杂的优化问题进行分解（decomposition）。谈到整数规划问题的分解主流的就三种方法：
1 Benders decomposition (主要思想是行生成+割平面方法)
2 Dantzig-Wolfe decomposition (主要思想其实就是列生成)
3 Lagrangian decomposition (主要思想是 Lagrangian relaxation)
前两种方法今天不是我们的主角，下面我们主要介绍 Lagrangian relaxation

1 拉格朗日松弛求解整数规划问题基本思想回顾

1.1 把约束放到目标上去 达到化繁为简的目的
这个我们就不具体讲了，因为基本上教科书上都会讲到。
1.2 尽量松弛掉linked/coupling constraints 让不可分问题变为可分问题
有如下优化问题
$\underset{x,y}{\min }{c}^{T}x+d^{T}y$
$A_{1}x=b_1$ (A1)
$A_{2}y=b_2$ (A2)
$A_{3}x+A_{4}y=b_3$ (A3)
决策变量有两大类 $x,y\in D$，约束（A1）和约束（A2）分别只和 $x,y$有关，而约束（A3）是一个linked/coupling constraints，它把$x,y$两类变量耦合在了一起。
若采用Lagrangian relaxation 把约束（A3）松弛到目标函数得到松弛问题：
$\underset{x,y}{\min }{c}^{T}x+d^{T}y+{\lambda }^T(A_{3}x+A_{4}y?b_3)$
$A_{1}x=b_1$ (B1)
$A_{2}y=b_2$ (B2)
此时由于linked/coupling constraints被处理掉了，因此在松弛问题中决策变量 $x,y$实现了解耦可以分为两个子问题分别对$x,y$求优化
子问题1：
$\underset{x}{\min }{c}^{T}x+{\lambda }^T{A}_{3}x$
$A_{1}x=b_1$
子问题2：
$\underset{y}{\min }{d}^{T}y+{\lambda }^T{A}_{4}y$
$A_{2}y=b_2$
为了后边解释方便我这边直接列出其对偶问题：
对偶函数：$q\left(\lambda \right)=\underset{A_{1}x=b_1,A_{2}y=b_2}{\min }{c}^{T}x+d^{T}y+{\lambda }^T\left(A_{3}x+A_{4}y?b_3\right)$
对偶问题：$\underset{\lambda }{\max }q\left(\lambda \right)$

2 为什么要用拉格朗日松弛？

1 拉格朗日松弛之后得到的问题下界要大于等于 直接将0,1整数变量松弛为[0,1]连续变量的下界。一言以蔽之就是 拉格朗松弛得到的解的质量要好。
2 采用拉格朗日松弛之后我们会得到拉格朗日乘子，拉格朗日乘子反应了对偶的信息，那么用这个这个对偶信息可以做很多很多事情（此处省略一万字）。
3 如果问题具有linked/coupling constraints那么采用拉格朗日松弛之后可以使得问题被分解，进而大大简化原问题。

3 拉格朗日迭代算法及收敛性分析

拉格朗日松弛求解整数规划问题的主要包括两大模块，1是松弛后的子问题求解；2是拉格朗日乘子的更新（对偶问题求解）
1松弛后的子问题求解
松弛后的子问题求解，一般要根据子问题的具体性质来确定求解方式。在整数规划问题中采用拉格朗日松弛后的子问题一般最好不要是NP-hard的或者一般问题规模很小可以进行暴力求解，否则子问题的求解还需要进一步特殊处理。所以一般我们求解子问题可以直接让求解器（Gurobi\Cplex等）直接求解即可，也可以根据具体情况采用动态规划，线性规划等手段帮助求解。
2拉格朗日乘子更新
拉格朗日乘子更新本质上是要求 对偶问题。由于对偶问题是一个非光滑凸优化问题，我们这里采用次梯度算法来更新拉格朗日乘子。如下所示：
${\lambda }_{k+1}={\lambda }_k+{\alpha }_k{g}_k$
其中$\lambda$ 为乘子，${\alpha }_k$,$g_k$分别为第k步迭代的步长和次梯度方向。其步长选取公式为
$0<{\alpha }_k<\frac{2\left(q^{?}?q\left({\lambda }_k\right)\right)}{\Vert g_k{\Vert }^2}$
次梯度选取公式为（次梯度的选择不唯一，这里只给出常见的一种）
$g_k=A_3{x}_k+A_4{y}_k?b_3$
松弛后的子问题求解和拉格朗日乘子更新分别反复迭代进行下去就构成了算法的流程。如下图所示：

只要按照公式（1-3）的方式选择步长和乘子迭代方向就可以保证收敛性，相关收敛性证明过程可以参考文献[3]的476页。

4 拉格朗日松弛迭代算法需要注意的问题

4.1 次梯度算法步长公式有大学问
公式（2）给出的步长选取范围
$0<{\alpha }_k<\frac{2\left(q^{?}?q\left({\lambda }_k\right)\right)}{\Vert g_k{\Vert }^2}$
其中 $q^{?}$是对偶问题最优解，实际上$q^{?}$是一个未知的量。我们整个算法就是为了求$q^{?}$那在求解过程中怎么会知道$q^{?}$?所以这个步长范围的公式看似美丽实际上不太实用，同时步长的选择其实对最终的求解效果有着举足轻重的影响。

目前已经有很多很多的paper着力于改进这个问题。一个思路是想办法估计$q^{?}$。另一个思路就是Surrogate Lagrangian Relaxation，具体可以参考文献[1].

4.2 尽量少松弛约束
平常在讲解拉格朗日松弛理论的时候是习惯把所有约束放到目标函数上去，以方便推导和讲解的过程，但是在实操过程中，松弛约束条件的做法实际上是没有办法的办法。能不松弛还是尽量还是不要松弛，能少松弛就少松弛一点。千万不要看着教科书上都是给全松弛了就一股脑全松弛掉约束条件，那样肯定不会得到高质量的解的。
4.3 松弛后的子问题要尽量简单
就像之前提到的拉格朗日松弛后的子问题一般最好不要是NP-hard的或者是问题规模很小即使是NP-hard也可以暴力求解的。
所以如果松弛后的子问题还是规模比较大的NP-hard问题，那么就需要一些特殊方法来进一步处理，那么有可能这个优化问题难度是比较大的，采用Lagrangian Relaxation未必能得到好的效果。
4.4 松弛后的问题尽量保证具备可分性
面对整数规划问题的求解很重要的一个思想就是通过decomposition来降低问题规模。如果所求解的优化问题没有 linked/coupling constraints 这种特性，可以在松弛linked/coupling constraints后将优化问题变成可分的性质的话，则该优化问题不推荐直接使用拉格朗日松弛来求解。
4.5 拉格朗日松弛得到的是一个最优解的下界（很可能是不可行解）
值得注意的是拉格朗日松弛得到的同样是一个原问题的下界，即不可行解。所以在完成拉格朗日松弛算法后我们还需要一个启发式算法把得到的这个不可行解修复为可行解。

5 拉格朗日松弛算法求解应用实例

要求解的简单二次整数规划问题实例，本例子来自于参考文献[1]的第一个实验
$\min 0.5x_1^2+0.1x_2^2+0.5x_3^2+0.1x_4^2+0.5x_5^2+0.1x_6^2$
$s.t.48?x_1+0.2x_2?x_3+0.2x_4?x_5+0.2x_6\le 0$
$250?5x_1+x_2?5x_3+x_4?5x_5+x_6\le 0$
x 1 , x 2 , x 3 , x 4 , x 5 , x 6 ∈ Z + ∪ { 0 } {x_1,x_2,x_3,x_4,x_5,x_6}\in Z_+\cup\left\{ 0 \right\} x1?,x2?,x3?,x4?,x5?,x6?∈Z+?∪{0}

将两个约束松弛到目标函数上来可得
$L(x_1,x_2,x_3,x_4,x_5,x_6,{\lambda }_1,{\lambda }_2)$
$=0.5x_1^2+0.1x_2^2+0.5x_3^2+0.1x_4^2+0.5x_5^2+0.1x_6^2{x}_5+0.2x_6$
$+{\lambda }_1(48?x_1+0.2x_2?x_3+0.2x_4?x_5+0.2x_6)$
$+{\lambda }_2(250?5x_1+x_2?5x_3+x_4?5x_5+x_6)$

代码主流程如下所述：

max_itertimes = 30
for i in range(max_itertimes):
    subproblem.compute_obj(dualproblem.lamd)
    subproblem.solve()
    subproblem.report(i)
    dualproblem.compute_subgradients(subproblem)
    dualproblem.compute_stepsize(subproblem)
    dualproblem.update_lamd()
    dualproblem.report(subproblem)
print("optimal solution = ", subproblem.opt_solution)
print("subgradients = ", dualproblem.subgradients)

subproblem为松弛问题子问题的对象，dualproblem为对偶问题的对象，然后我们先求解松弛问题 subproblem.solve()，然后再求解对偶问题，求解对偶包括 计算次梯度（compute_subgradients），计算步长(compute_stepsize)和更新乘子(update_lamd)
本案例完整代码参考：https://github.com/WenYuZhi/lagrangianRelaxationQIP

我们（一作我师妹，三作是我，我主要负责了代码开发）最近也做了一个关于Surrogate Lagrangian relaxation求解生产调度的问题见参考文献[2]，后期我们也会把代码开源供大家学习参考。

参考文献：
【1】Bragin M A, Luh P B, Yan J H, et al. Convergence of the Surrogate Lagrangian Relaxation Method[J]. Journal of Optimization Theory and Applications, 2015, 164(1): 173-201.
【2】Cui H, Luo X, Wang Y, et al. Scheduling of steelmaking-continuous casting process using deflected surrogate Lagrangian relaxation approach and DC algorithm[J]. Computers & Industrial Engineering, 2020.
【3】Bertsekas D P. 非线性规划[J]. 清华大学出版社，2013.