导读

拉格朗日松弛（Lagrangian Relaxation）是Held和Krap在研究旅行商问题时提出的一种优化方法。算法的基本原理是通过引入拉格朗日松弛算子将约束条件中造成问题“求解困难”的约束吸收到目标函数中, 并保持目标函数的线性形式，从而使松弛后的问题容易求解。本文将介绍拉格朗日松弛算法的基本理论及算法流程，并以约束最短路问题为例使用Python语言进行编程求解。具体内容包括：

1.介绍拉格朗日松弛算法的基本理论和算法流程

2.介绍基于拉格朗日松弛算法的约束最短路问题

3.介绍约束最短路问题的Python代码

拉格朗日松弛算法

1.1 拉格朗日函数及拉格朗日对偶

考虑以下极小化问题：

若由于g_i(x)≥b_i,i=1,2,......,m的存在使得问题求解更加困难，则称其为困难约束。而xX称为简单约束。为便于求解，可以通过引入拉格朗日乘子作为困难约束的惩罚项将困难约束写到目标函数中：

上式也可以写作

我们也称之为拉格朗日函数。在拉格朗日函数的基础上，我们可以得到拉格朗日对偶：

根据弱对偶定理，每组拉格朗日乘子（即拉格朗日对偶问题的一个解）所对应的目标函数值即为原始问题最优解的一个下界。如果获得拉格朗日对偶问题更高质量的解（即更佳的拉格朗日乘子），就可以使更加接近原始问题的最优解。因此拉格朗日松弛算法的基本思路就是通过更新拉格朗日乘子，从而获得原问题更好地下界。

1.2 拉格朗日乘子的确定（次梯度法迭代法）

为了能够获得最佳的一组拉格朗日乘子，我们采用次梯度（Subgradient）迭代法，即通过更新拉格朗日乘子，逐渐提升下界。次梯度算法如下：

Step0：任选初始拉格朗日乘子，t=1；

Step1：对，求解，得到其最优解；

Step2：利用计算次梯度；

Step3：若满足某种停止原则（下文列出常用的几种），则算法停止，获得原始问题的下界，否则，转Step4；

Step4：令t=t+1，更新乘子，其中，转Step1。

1.3 停止原则

①迭代次数不超过T。这是一种最为简单的原则，但解的质量无法保证。

②=0，这是最为理想的状态，此时达到拉格朗日对偶问题的最优解。在实际计算中，由于问题的复杂性和机器精度误差，这样的结果很难达到，常用（为给定的非负数）来代替。

③。我们在使用拉格朗日松弛算法的迭代过程中，可以通过更新拉格朗日乘子逐步得到问题更好的下界。而上界的确定往往要借助于启发式算法，即通过启发式算法将松弛问题所得的解转换为可行解所对应的目标函数值作为问题的上界UB。当上界UB和当前下界足够接近时，即当时，表示已经得到原问题的（近似）最优解，可以终止算法，并证明启发式算法所获得的可行解具有较高的质量。此外，在设计获取问题可行解的算法时也可利用拉格朗日松弛算法迭代过程中所获得的信息，从而在拉格朗日算法每次迭代后，能够获得一个更高质量（即更小）的上界。

④或目标值在规定步数内的变化不超过事先给定的值，则认为目标值不可能再变化，停止运算。

约束最短路问题

2.1 问题描述

考虑如图1所示网络，其中每条弧上标注了相应的花费和行程时间。期望找到由源节点1到汇聚节点6的最短路径，且保证行程时间不能超过T（设置为15）。

2.2 模型构建

根据2.1问题描述，可以构建模型如下：

其中决策变量=1表示弧(i,j)被选择，否则=0。

2.3 拉格朗日松弛模型

将模型中的时间约束借助拉格兰朗日乘子放入目标中，可以将问题转化为最短路问题，松弛后的模型为：

2.4 迭代过程参数设置

设置初始拉格朗日乘子，在计算迭代步长时设置，UB=20，设置停止准则为次梯度=0或者迭代次数t=500。

Python代码实现

3.1 gurobi包的导入及参数设置

导入gurobi包，利用gurobi来构建模型并求解。

from gurobipy import *

输入构建模型所需数据，设置控制迭代过程的参数。

Arc, Cost = multidict({(1,2):1, #费用

        (1,3):10,

        (2,4):1,

        (2,5):2,

        (3,2):1,

        (3,4):5,

        (3,5):12,

        (4,5):10,

        (4,6):1,

        (5,6):2,})



Time = {(1,2):10,               #时间

        (1,3):3,

        (2,4):1,

        (2,5):3,

        (3,2):2,

        (3,4):7,

        (3,5):3,

        (4,5):1,

        (4,6):7,

        (5,6):2,}

n = 6           #节点个数

T = 15          #约束最短路的时间限制

lam = 0         #初始拉格朗日乘子

UB = 20         #上界

beta = 0.1      #迭代步长的倍数

3.2 拉格朗日松弛模型构建

利用gurobi包来构建拉格朗日松弛模型并进行求解。

mm = Model("约束最短路")

x = mm.addVars(Arc, vtype = GRB.BINARY)

mm.setObjective(quicksum(Cost*x for i in Arc) + lam*(quicksum(Time*x for i in Arc)-T))

for i in range(1,n+1):

    if i == 1:

        mm.addConstr(quicksum(x for i in Arc.select(i,'*')) == 1)

    elif i == n:

        mm.addConstr(quicksum(x for i in Arc.select('*',i)) == 1)

    else:

        mm.addConstr(quicksum(x for i in Arc.select(i,'*'))-quicksum(x for i in Arc.select('*',i)) == 0)

mm.optimize()

3.3 使用次梯度法进行迭代

获取当前求得的下界值，并计算当前的次梯度。

DB = mm.objVal     #获取当前下界值  

subgradient = quicksum(Time*x.x for i in Arc) - T  #计算次梯度

判断算法是否停止。

if subgradient.getValue() == 0:     #若次梯度为0，则停止

    break

if t >= 500:                        #若迭代次数达到500次，则停止

    break

计算步长，更新。

theta = beta * (UB - DB) / subgradient.getValue()**2 #计算步长

lam = max(0, lam + theta * subgradient.getValue())  #更新拉格朗日乘子

在更新后重新求解拉格朗日松弛模型，并使迭代次数加1。

mm.setObjective(quicksum(Cost*x for i in Arc) + lam*(quicksum(Time*x for i in Arc)-T))

mm.optimize()

t=t+1

3.4 计算结果

迭代过程如下表所示：

获得最优解为(1,2,5,6)，最小费用为5，总花费时间为15，满足时间要求。注意，当改变时间限制的取值时，可能无法直接得到原问题的最优解，而只能收敛到一个比较好的下界值。此时可以通过其他的方法（比如启发式算法）来得到可行解，并根据拉格朗日松弛得到的下界值来判断可行解的质量。

参考文献

[1] Held, M., & Karp, R. M. (1970). The traveling-salesman problem and minimum spanning trees. Operations Research, 18(6), 1138-1162.

[2] Held, M., & Karp, R. M. (1971). The traveling-salesman problem and minimum spanning trees: Part II. Mathematical programming, 1(1), 6-25.

[3] Bazaraa, M. S., Sherali, H. D., & Shetty, C. M. (2013). Nonlinear programming: theory and algorithms. John Wiley & Sons.

---

本文转载自微信公众号：【交通与优化】（ID:OR_Transportation）

作者：王鸿杨

导读

拉格朗日松弛（Lagrangian Relaxation）是Held和Krap在研究旅行商问题时提出的一种优化方法。算法的基本原理是通过引入拉格朗日松弛算子将约束条件中造成问题“求解困难”的约束吸收到目标函数中, 并保持目标函数的线性形式，从而使松弛后的问题容易求解。本文将介绍拉格朗日松弛算法的基本理论及算法流程，并以约束最短路问题为例使用Python语言进行编程求解。具体内容包括：

1.介绍拉格朗日松弛算法的基本理论和算法流程

2.介绍基于拉格朗日松弛算法的约束最短路问题

3.介绍约束最短路问题的Python代码

拉格朗日松弛算法

1.1 拉格朗日函数及拉格朗日对偶

考虑以下极小化问题：

若由于g_i(x)≥b_i,i=1,2,......,m的存在使得问题求解更加困难，则称其为困难约束。而xX称为简单约束。为便于求解，可以通过引入拉格朗日乘子作为困难约束的惩罚项将困难约束写到目标函数中：

上式也可以写作

我们也称之为拉格朗日函数。在拉格朗日函数的基础上，我们可以得到拉格朗日对偶：

根据弱对偶定理，每组拉格朗日乘子（即拉格朗日对偶问题的一个解）所对应的目标函数值即为原始问题最优解的一个下界。如果获得拉格朗日对偶问题更高质量的解（即更佳的拉格朗日乘子），就可以使更加接近原始问题的最优解。因此拉格朗日松弛算法的基本思路就是通过更新拉格朗日乘子，从而获得原问题更好地下界。

1.2 拉格朗日乘子的确定（次梯度法迭代法）

为了能够获得最佳的一组拉格朗日乘子，我们采用次梯度（Subgradient）迭代法，即通过更新拉格朗日乘子，逐渐提升下界。次梯度算法如下：

Step0：任选初始拉格朗日乘子，t=1；

Step1：对，求解，得到其最优解；

Step2：利用计算次梯度；

Step3：若满足某种停止原则（下文列出常用的几种），则算法停止，获得原始问题的下界，否则，转Step4；

Step4：令t=t+1，更新乘子，其中，转Step1。

1.3 停止原则

①迭代次数不超过T。这是一种最为简单的原则，但解的质量无法保证。

②=0，这是最为理想的状态，此时达到拉格朗日对偶问题的最优解。在实际计算中，由于问题的复杂性和机器精度误差，这样的结果很难达到，常用（为给定的非负数）来代替。

③。我们在使用拉格朗日松弛算法的迭代过程中，可以通过更新拉格朗日乘子逐步得到问题更好的下界。而上界的确定往往要借助于启发式算法，即通过启发式算法将松弛问题所得的解转换为可行解所对应的目标函数值作为问题的上界UB。当上界UB和当前下界足够接近时，即当时，表示已经得到原问题的（近似）最优解，可以终止算法，并证明启发式算法所获得的可行解具有较高的质量。此外，在设计获取问题可行解的算法时也可利用拉格朗日松弛算法迭代过程中所获得的信息，从而在拉格朗日算法每次迭代后，能够获得一个更高质量（即更小）的上界。

④或目标值在规定步数内的变化不超过事先给定的值，则认为目标值不可能再变化，停止运算。

约束最短路问题

2.1 问题描述

考虑如图1所示网络，其中每条弧上标注了相应的花费和行程时间。期望找到由源节点1到汇聚节点6的最短路径，且保证行程时间不能超过T（设置为15）。

2.2 模型构建

根据2.1问题描述，可以构建模型如下：

其中决策变量=1表示弧(i,j)被选择，否则=0。

2.3 拉格朗日松弛模型

将模型中的时间约束借助拉格兰朗日乘子放入目标中，可以将问题转化为最短路问题，松弛后的模型为：

2.4 迭代过程参数设置

设置初始拉格朗日乘子，在计算迭代步长时设置，UB=20，设置停止准则为次梯度=0或者迭代次数t=500。

Python代码实现

3.1 gurobi包的导入及参数设置

导入gurobi包，利用gurobi来构建模型并求解。

from gurobipy import *

输入构建模型所需数据，设置控制迭代过程的参数。

Arc, Cost = multidict({(1,2):1, #费用

        (1,3):10,

        (2,4):1,

        (2,5):2,

        (3,2):1,

        (3,4):5,

        (3,5):12,

        (4,5):10,

        (4,6):1,

        (5,6):2,})



Time = {(1,2):10,               #时间

        (1,3):3,

        (2,4):1,

        (2,5):3,

        (3,2):2,

        (3,4):7,

        (3,5):3,

        (4,5):1,

        (4,6):7,

        (5,6):2,}

n = 6           #节点个数

T = 15          #约束最短路的时间限制

lam = 0         #初始拉格朗日乘子

UB = 20         #上界

beta = 0.1      #迭代步长的倍数

3.2 拉格朗日松弛模型构建

利用gurobi包来构建拉格朗日松弛模型并进行求解。

mm = Model("约束最短路")

x = mm.addVars(Arc, vtype = GRB.BINARY)

mm.setObjective(quicksum(Cost*x for i in Arc) + lam*(quicksum(Time*x for i in Arc)-T))

for i in range(1,n+1):

    if i == 1:

        mm.addConstr(quicksum(x for i in Arc.select(i,'*')) == 1)

    elif i == n:

        mm.addConstr(quicksum(x for i in Arc.select('*',i)) == 1)

    else:

        mm.addConstr(quicksum(x for i in Arc.select(i,'*'))-quicksum(x for i in Arc.select('*',i)) == 0)

mm.optimize()

3.3 使用次梯度法进行迭代

获取当前求得的下界值，并计算当前的次梯度。

DB = mm.objVal     #获取当前下界值  

subgradient = quicksum(Time*x.x for i in Arc) - T  #计算次梯度

判断算法是否停止。

if subgradient.getValue() == 0:     #若次梯度为0，则停止

    break

if t >= 500:                        #若迭代次数达到500次，则停止

    break

计算步长，更新。

theta = beta * (UB - DB) / subgradient.getValue()**2 #计算步长

lam = max(0, lam + theta * subgradient.getValue())  #更新拉格朗日乘子

在更新后重新求解拉格朗日松弛模型，并使迭代次数加1。

mm.setObjective(quicksum(Cost*x for i in Arc) + lam*(quicksum(Time*x for i in Arc)-T))

mm.optimize()

t=t+1

3.4 计算结果

迭代过程如下表所示：

获得最优解为(1,2,5,6)，最小费用为5，总花费时间为15，满足时间要求。注意，当改变时间限制的取值时，可能无法直接得到原问题的最优解，而只能收敛到一个比较好的下界值。此时可以通过其他的方法（比如启发式算法）来得到可行解，并根据拉格朗日松弛得到的下界值来判断可行解的质量。

参考文献

[1] Held, M., & Karp, R. M. (1970). The traveling-salesman problem and minimum spanning trees. Operations Research, 18(6), 1138-1162.

[2] Held, M., & Karp, R. M. (1971). The traveling-salesman problem and minimum spanning trees: Part II. Mathematical programming, 1(1), 6-25.

[3] Bazaraa, M. S., Sherali, H. D., & Shetty, C. M. (2013). Nonlinear programming: theory and algorithms. John Wiley & Sons.

---

本文转载自微信公众号：【交通与优化】（ID:OR_Transportation）

作者：王鸿杨