大规模线性规划的对偶问题
对偶理论简介
SPP的对偶
MCNF的对偶
下面是公众号的稿子
运筹学修炼日记：如何优雅地写出大规模线性规划的对偶
对偶理论Duality Theory在运筹学数学规划部分占据着举足轻重的地位，也属于比较高阶的理论。Duality Theory在精确算法设计中也经常用到，在Robust Optimization等涉及到多层规划Multilevel的问题中，也有非常广泛的应用，很多时候可以化腐朽为神奇。尤其在Robust Optimization中，有些问题可以巧妙的将内层inner level的模型转化成LP，从而可以通过对偶，将双层bi-level的模型，转化成单阶段single level的模型，从而用单层的相关算法来求解Robust Optimization问题。
今天我们就来看看，在实际的科研当中，遇到的一些稍微复杂一点的LP，我们如何写出其对偶问题。
实际上在一些顶刊中，例如transportation Science等，比较近期的文章，也时不时会看到这样的操作。这个操作其实并不是抬手就能搞定的，很多时候需要反复修改，才能将对偶问题正确的写出来。
据我所知，我似乎是第一个写这样博文的博主。(如果有比我更早的，请告知我掐了这段)
先来看一个比较容易的线性规划问题：
max
Z
=
2
x
1
+
3
x
2
5
x
1
+
4
x
2
⩽
170
→
y
1
2
x
1
+
3
x
2
⩽
100
→
y
2
x
1
,
x
1
⩾
0
\begin{aligned} \max \quad Z=&2x_1+3x_2 \\ &5x_1+4x_2\leqslant 170 \quad \rightarrow \quad \color{red} y_1 \\ &2x_1+3x_2\leqslant 100 \quad \rightarrow \quad \color{red} y_2 \\ &x_1,x_1\geqslant 0 \end{aligned}
maxZ=2x1+3x25x1+4x2⩽170→y12x1+3x2⩽100→y2x1,x1⩾0
其对偶问题比较容易写出：
min
W
=
170
y
1
+
100
y
2
5
y
1
+
2
y
2
⩾
2
4
y
1
+
3
y
2
⩾
3
y
1
,
y
1
⩾
0
\begin{aligned} \min \quad W=&170y_1+100y_2 \\ &5y_1+2y_2\geqslant 2 \\ &4y_1+3y_2\geqslant 3 \\ &y_1,y_1\geqslant 0 \end{aligned}
minW=170y1+100y25y1+2y2⩾24y1+3y2⩾3y1,y1⩾0
基本原则如图：
但是假如是最短路问题：
最短路问题
max ∑ e ∈ A d e x e ∑ e ∈ o u t ( i ) x e − ∑ e ∈ i n ( i ) x e = { 1 , if i = s − 1 , if i = t 0 , else 0 ⩽ x e ⩽ 1 , ∀ e ∈ A \begin{aligned} \max \quad & \sum_{e\in A}{d_ex_e} \\ &\sum_{e\in \mathrm{out}\left( i \right)}{x_e}-\sum_{e\in \mathrm{in}\left( i \right)}{x_e}=\begin{cases} 1,& \text{if}\,\,i=s\\ -1,& \text{if}\,\,i=t\\ 0,& \text{else}\\ \end{cases} \\ & 0 \leqslant x_e \leqslant 1 ,\qquad \forall e\in A \end{aligned} maxe∈A∑dexee∈out(i)∑xe−e∈in(i)∑xe=⎩⎪⎨⎪⎧1,−1,0,ifi=sifi=telse0⩽xe⩽1,∀e∈A
这里大括号里有几个条件判断，就不是那么容易了。也许这个还比较容易，那再看看这个。多商品流问题Multicommodity Network Flow Problem
多商品流问题Multicommodity Network Flow Problem
- K K K origin-destination pairs of nodes, ( s 1 , t 1 , d 1 ) , ( s 2 , t 2 , d 2 ) , ⋯ , ( s k , t k , d k ) (s_1, t_1, d_1), (s_2, t_2, d_2), \cdots, (s_k, t_k, d_k) (s1,t1,d1),(s2,t2,d2),⋯,(sk,tk,dk).
- d k d_k dk : demand, amount of flow that must be sent from s k s_k sk to t k t_k tk.
- u i j u_{ij} uij : capacity on ( i , j ) (i,j) (i,j) shared by all commodities
- c i j k c_{ij}^k cijk : cost of sending 1 unit of commodity k k k in ( i , j ) (i,j) (i,j)
- x i j k x_{ij}^k xijk : flow of commodity k k k in ( i , j ) (i,j) (i,j)
模型如下
min
∑
(
i
,
j
)
∈
A
∑
k
c
i
j
k
x
i
j
k
∑
j
x
i
j
k
−
∑
j
x
j
i
k
=
{
d
k
,
i
f
i
=
s
k
−
d
k
,
i
f
i
=
t
k
0
,
o
t
h
e
r
w
i
s
e
∑
k
x
i
j
k
⩽
u
i
j
,
∀
(
i
,
j
)
∈
A
x
i
j
k
⩾
0
,
∀
(
i
,
j
)
∈
A
,
k
∈
K
\begin{aligned} \min \quad &\sum_{\left( i,j \right) \in A}{\sum_k{c_{ij}^{k}x_{ij}^{k}}} \\ &\sum_j{x_{ij}^{k}}-\sum_j{x_{ji}^{k}}=\begin{cases} d_k,& \qquad \mathrm{if}\,\,\,\, i=s_k\\ -d_k,& \qquad \mathrm{if}\,\,\,\, i=t_k\\ 0,& \qquad \mathrm{otherwise}\\ \end{cases} \\ &\sum_k{x_{ij}^{k}}\leqslant u_{ij}, \qquad \forall \left( i,j \right) \in A \\ &x_{ij}^{k}\geqslant 0, \qquad \forall \left( i,j \right) \in A, k\in K \end{aligned}
min(i,j)∈A∑k∑cijkxijkj∑xijk−j∑xjik=⎩⎪⎨⎪⎧dk,−dk,0,ifi=skifi=tkotherwisek∑xijk⩽uij,∀(i,j)∈Axijk⩾0,∀(i,j)∈A,k∈K
现在呢？还能很容易的写出来吗？若果还能，大神请受小弟一拜！哈哈哈
能轻松写出来的非人类大神可以前走左拐去刷剧了，咱这些普通人就接着往下看把。
注意，上面的Multicommodity Flow Problem和Shortest Path Problem都是Linear Programming,可以对偶的。但是对于Integer Programming和Mixed Integer Programming来讲，是不能对偶的。这一点一定要搞清楚。
对于这种稍微复杂一些的LP，我们怎么能写出对偶还保证正确，可debug找错的呢？我的方法就是借助Excel+具体小算例。
借助Excel和具体小算例写出大规模LP的对偶
为了大家理解方便，我们不要直接去硬钢Multicommodity Network Flow(理论上搞定了Multicommodity Network Flow，其实就具备搞定大多数可以对偶的LP的潜力了).我们先以SPP来开个胃。
Dual Problem :Shortest Path Problem(最短路问题)
小算例
我们先来引入一个小算例。该算例来自参考文献[^2],我做了点修改。为了显示我比较认真，我还专门无聊用LaTeX+Tikz重新画了一个小花图：
(咋样，看着还舒服吧)
我们再来看一下SPP的模型：
max
∑
e
∈
A
d
e
x
e
∑
e
∈
o
u
t
(
i
)
x
e
−
∑
e
∈
i
n
(
i
)
x
e
=
{
1
,
if
i
=
s
−
1
,
if
i
=
t
0
,
else
0
⩽
x
e
⩽
1
,
∀
e
∈
A
\begin{aligned} \max \quad & \sum_{e\in A}{d_ex_e} \\ &\sum_{e\in \mathrm{out}\left( i \right)}{x_e}-\sum_{e\in \mathrm{in}\left( i \right)}{x_e}=\begin{cases} 1,& \text{if}\,\,i=s\\ -1,& \text{if}\,\,i=t\\ 0,& \text{else}\\ \end{cases} \\ & 0 \leqslant x_e \leqslant 1 ,\qquad \forall e\in A \end{aligned}
maxe∈A∑dexee∈out(i)∑xe−e∈in(i)∑xe=⎩⎪⎨⎪⎧1,−1,0,ifi=sifi=telse0⩽xe⩽1,∀e∈A
按照这个模型，我们手动把这个模型具体的写出来。为了之后的操作，我们直接写到Excel里。
Excel+小算例写出SPP的对偶问题
SPP模型如下：
我们把这个表叫Primal tabular其中，每一列代表一个变量
x
i
j
,
∀
(
i
,
j
)
∈
A
x_{ij}, \forall (i,j)\in A
xij,∀(i,j)∈A
- 第一行代表该条 ( i , j ) (i,j) (i,j)的距离
- 第二行代表变量 x i j , ∀ ( i , j ) ∈ A x_{ij}, \forall (i,j)\in A xij,∀(i,j)∈A
- 第一行和第二行就组成了目标函数 ∑ e ∈ A d e x e \sum_{e\in A}{d_ex_e} ∑e∈Adexe
- 第3-9行代表每个结点 ∀ i ∈ V \forall i \in V ∀i∈V的约束
- 最后一列代表每个约束的
Dual variable
OK，我们按照对偶的方法，将Primal tabular的RHS和Dual variabe拷贝，转置成2行，放在一个新表格(我们叫做Dual tabular)的头两行，然后将Primal tabular的整个约束系数矩阵拷贝，转置到Dual tabular头两行下面。再把原问题Primal tabular的distance行和min行拷贝，转置，放在Dual tabular的右面。
再把Dual tabular中改成max。为了明确Dual Problem中各个变量的符号（正负性）以及每个约束的符号(relation)，我们在Dual tabular中加入一行（就是第三行）,表示变量的符号。同时在Dual tabular约束矩阵后加入一列，表示约束的符号。
操作完就是这样的
按照上面那个关系图中的信息，我们可以确定，对偶变量
π
i
\pi_i
πi都是无约束的，我们用=表示，Dual Problem中的约束都是
⩽
\leqslant
⩽的。这样，对偶就完成了。
但是，这还是一个具体的算例的Dual，我们需要将这个具体的算例，通过提取信息整理，化成一个general的公式形式。
将Excel中的Dual tabular转化成公式形式
我们观察上图，每一行都对应一条弧
(
i
,
j
)
∈
A
(i,j)\in A
(i,j)∈A,例如第一行是
(
1
,
2
)
(1, 2)
(1,2)，第二行是
(
1
,
4
)
(1,4)
(1,4)等。可以看到，对应出发点的变量系数全是1，对应终点的系数全是-1，无一例外，因此，我们可以断定，这个约束可以这么写：
π
i
−
π
j
⩽
d
i
j
,
∀
(
i
,
j
)
∈
A
\pi _i-\pi _j \leqslant d_{ij}, \qquad \forall \left( i,j \right) \in A
πi−πj⩽dij,∀(i,j)∈A
结合目标函数，以及变量的符号，我们可以写出SPP的对偶问题：
max
π
s
−
π
t
π
i
−
π
j
⩽
d
i
j
,
∀
(
i
,
j
)
∈
A
π
i
free
\begin{aligned} \max \qquad &\pi _s-\pi _t \\ &\pi _i-\pi _j \leqslant d_{ij}, \qquad \forall \left( i,j \right) \in A \\ &\pi _i\,\,\,\,\text{free} \end{aligned}
maxπs−πtπi−πj⩽dij,∀(i,j)∈Aπifree
大功告成，怎么样，有没有点内味了？
接下来，我们啃一个稍微难啃一些的骨头Multicommodity Network Flow Problem.
Dual Problem :Multicommodity Network Flow Problem(最短路问题)
这个问题相比SPP难度还是大挺多的。我们首先上数学模型。
min ∑ ( i , j ) ∈ A ∑ k c i j k x i j k ∑ j x i j k − ∑ j x j i k = { d k , i f i = s k − d k , i f i = t k 0 , o t h e r w i s e ∑ k x i j k ⩽ u i j , ∀ ( i , j ) ∈ A x i j k ⩾ 0 , ∀ ( i , j ) ∈ A , k ∈ K \begin{aligned} \min \quad &\sum_{\left( i,j \right) \in A}{\sum_k{c_{ij}^{k}x_{ij}^{k}}} \\ &\sum_j{x_{ij}^{k}}-\sum_j{x_{ji}^{k}}=\begin{cases} d_k,& \qquad \mathrm{if}\,\,\,\, i=s_k\\ -d_k,& \qquad \mathrm{if}\,\,\,\, i=t_k\\ 0,& \qquad \mathrm{otherwise}\\ \end{cases} \\ &\sum_k{x_{ij}^{k}}\leqslant u_{ij}, \qquad \forall \left( i,j \right) \in A \\ &x_{ij}^{k}\geqslant 0, \qquad \forall \left( i,j \right) \in A, k\in K \end{aligned} min(i,j)∈A∑k∑cijkxijkj∑xijk−j∑xjik=⎩⎪⎨⎪⎧dk,−dk,0,ifi=skifi=tkotherwisek∑xijk⩽uij,∀(i,j)∈Axijk⩾0,∀(i,j)∈A,k∈K
看看吧，又有什么
i
f
i
=
s
k
\mathrm{if}\,\, i=s_k
ifi=sk之类的，变量还是
x
i
j
k
x_{ij}^k
xijk,你尝试自己先写一下，是不是觉得脑瓜子嗡嗡的，哈哈哈。
具体算例
都不是事儿，咱一起来刚一下。同样的把文献[^2]中的算例原模原样搬过来看看（当然图还是我自己画的）
Excel+小算例写出MNF的原模型
我们考虑有两个commodity:
commodity = [[1, 7, 25], # s_i, d_i, demand
[2, 6, 2]
]
本来想把代码也放上的，感觉太多了，有需要的话，大家私信我，我在修改把代码放上来。
然后我们按照模型和算例网络结构，把模型具体的写出来，如下图所示
第二行的变量x_1_2_0就代表 x i j k x_{ij}^k xijk，其中0代表 k k k。
这个看上去不太有规律，我们按照commodity k k k把上面的表格整一下，变成：
现在看上去就比较清楚了。我们仍然把这个表格叫做Primal Tabular.
接下来我们按照同样的方法，根据Primal Tabular生成Dual Tabular,如下图
将Excel中的Dual tabular转化成公式形式
为了区分 u u u和 μ \mu μ，表格中的mu我就用 λ \lambda λ代替了，因为表格中写mu省地儿。
所以大家注意 λ i j \lambda _{ij} λij就是上面表格中的mu
max ∑ k ∈ K d k ( π i = s k k − π i = t k k ) + ∑ ( i , j ) ∈ A u i j λ i j π i k − π j k + λ i j ⩽ c i j k , ∀ k ∈ K , ∀ ( i , j ) ∈ A π i k free , ∀ k ∈ K , ∀ i ∈ V λ i j ⩽ 0 , ∀ ( i , j ) ∈ A \begin{aligned} \max & \sum_{k\in K}{d_k\left( \pi _{i=s_k}^{k}-\pi _{i=t_k}^{k} \right)}+\sum_{\left( i,j \right) \in A}{u_{ij}\lambda _{ij}} \\ &\pi _{i}^{k}-\pi _{j}^{k}+\lambda _{ij}\leqslant c_{ij}^{k}, \qquad \forall k\in K,\forall \left( i,j \right) \in A \\ &\pi _{i}^{k}\,\,\,\,\text{free}, \qquad \forall k\in K,\forall i\in V \\ &\lambda _{ij}\leqslant 0, \qquad \forall \left( i,j \right) \in A \end{aligned} maxk∈K∑dk(πi=skk−πi=tkk)+(i,j)∈A∑uijλijπik−πjk+λij⩽cijk,∀k∈K,∀(i,j)∈Aπikfree,∀k∈K,∀i∈Vλij⩽0,∀(i,j)∈A
当然了，按照国际惯例（搞OR大佬的惯例），我们还是跟之前我写的讲SPP对偶的博文
https://blog.csdn.net/HsinglukLiu/article/details/107834197
中的操作一样：
- 将所有对偶变量 π i k \pi _{i}^{k} πik取相反数
- 把原约束中 π i k − π j k \pi _{i}^{k}-\pi _{j}^{k} πik−πjk改成 π j k − π i k \pi _{j}^{k}-\pi _{i}^{k} πjk−πik
- 将 π i = s k k \pi _{i=s_k}^{k} πi=skk设置成0，也就是 π i = s k k = 0 \pi _{i=s_k}^{k}=0 πi=skk=0
这三个隐含小动作，大佬是不会在论文里面写的，要是没仔细钻研，你一般会一头雾水。
OK，按照国际惯例操作完后，最终Multicommodity Network Flow Problem模型的Dual Problem就变成了下面的样子
max
∑
k
∈
K
d
k
π
i
=
t
k
k
+
∑
(
i
,
j
)
∈
A
u
i
j
λ
i
j
π
j
k
−
π
i
k
+
λ
i
j
⩽
c
i
j
k
,
∀
k
∈
K
,
∀
(
i
,
j
)
∈
A
π
i
k
free
,
π
s
k
k
=
0
,
∀
k
∈
K
,
∀
i
∈
V
λ
i
j
⩽
0
,
∀
(
i
,
j
)
∈
A
\begin{aligned} \max & \sum_{k\in K}{d_k\pi _{i=t_k}^{k}}+\sum_{\left( i,j \right) \in A}{u_{ij}\lambda _{ij}} \\ &\pi _{j}^{k}-\pi _{i}^{k}+\lambda _{ij}\leqslant c_{ij}^{k}, \qquad \forall k\in K,\forall \left( i,j \right) \in A \\ &\pi _{i}^{k}\,\,\,\,\text{free}, \pi _{s_k}^{k} = 0, \qquad \forall k\in K,\forall i\in V \\ &\lambda _{ij}\leqslant 0, \qquad \forall \left( i,j \right) \in A \end{aligned}
maxk∈K∑dkπi=tkk+(i,j)∈A∑uijλijπjk−πik+λij⩽cijk,∀k∈K,∀(i,j)∈Aπikfree,πskk=0,∀k∈K,∀i∈Vλij⩽0,∀(i,j)∈A
OK,所有的动作都完成了。一块硬骨头啃完了。
Python调用Gurobi求解Multicommodity Network Flow Problem （仅原问题）
最后再附上求解这个问题的Python代码（对偶问题的不想写了）
from gurobipy import *
import pandas as pd
import numpy as np
from pandas import Series, DataFrame
import math
Arcs = {'1,2': [15, 15] # cost flow
,'1,4': [25, 25]
,'1,3': [45, 45]
,'2,5': [30, 60]
,'2,4': [2, 2]
,'5,7': [2, 2]
,'4,7': [50, 100]
,'4,3': [2, 2]
,'3,6': [25, 50]
,'6,7': [1, 1]
}
Arcs
Nodes = [1, 2, 3, 4, 5, 6, 7]
commodity = [[1, 7, 25], # s_i, d_i, demand
[2, 6, 2]
]
model = Model('MultiCommodity')
# add variables
X = {}
for key in Arcs.keys():
for k in range(len(commodity)):
key_x = key + ',' + str(k)
X[key_x] = model.addVar(lb=0
,ub=Arcs[key][1]
,vtype=GRB.CONTINUOUS
,name= 'x_' + key_x
)
# add objective function
obj = LinExpr(0)
for key in Arcs.keys():
for k in range(len(commodity)):
key_x = key + ',' + str(k)
obj.addTerms(Arcs[key][0], X[key_x])
model.setObjective(obj, GRB.MINIMIZE)
# constraints 1
for k in range(len(commodity)):
for i in Nodes:
lhs = LinExpr(0)
for key_x in X.keys():
# nodes = key_x.split(',')
if(i == (int)(key_x.split(',')[0]) and k == (int)(key_x.split(',')[2])):
lhs.addTerms(1, X[key_x])
if(i == (int)(key_x.split(',')[1]) and k == (int)(key_x.split(',')[2])):
lhs.addTerms(-1, X[key_x])
if(i == commodity[k][0]):
model.addConstr(lhs == commodity[k][2], name='org_, ' + str(i) + '_' + str(k))
elif(i == commodity[k][1]):
model.addConstr(lhs == -commodity[k][2], name='des_, ' + str(i) + '_' + str(k))
else:
model.addConstr(lhs == 0, name='inter_, ' + str(i) + '_' + str(k))
# constraints 2
for key in Arcs.keys():
lhs = LinExpr(0)
for k in range(len(commodity)):
key_x = key + ',' + str(k)
lhs.addTerms(1, X[key_x])
model.addConstr(lhs <= Arcs[key][1], name = 'capacity_, ' + key)
model.write('Multicommodity_model.lp')
model.optimize()
for var in model.getVars():
if(var.x > 0):
print(var.varName, '\t', var.x)
dual = model.getAttr("Pi", model.getConstrs())
原问题求解结果如下：
Solved in 0 iterations and 0.01 seconds
Optimal objective 1.873000000e+03
x_1,2,0 2.0
x_1,4,0 22.0
x_1,3,0 1.0
x_2,5,0 2.0
x_2,4,1 2.0
x_5,7,0 2.0
x_4,7,0 22.0
x_4,3,1 2.0
x_3,6,0 1.0
x_3,6,1 2.0
x_6,7,0 1.0
对偶问题求解
后记
硕士的时候搞这个搞了几天，还请教了我师兄挺多。师兄的研究Robust Service Network Design的文章里也用到了类似这样问题的对偶，发了Transportation Science，我把文章也贴在这里，欢迎大家去读一读，做的非常好[^3]。可以看到，这样的技巧在科研中还是有用武之地的。
[1] :Garg, N., & Koenemann, J. (2007). Faster and simpler algorithms for multicommodity flow and other fractional packing problems. SIAM Journal on Computing, 37(2), 630-652https://doi.org/10.1137/S0097539704446232
[2]:Cappanera, P., & Scaparra, M. P. (2011). Optimal allocation of protective resources in shortest-path networks. Transportation Science, 45(1), 64-80.
http://dx.doi.org/10.1287/trsc.1100.0340
[3]:Wang, Z., & Qi, M. (2020). Robust service network design under demand uncertainty. Transportation Science, 54(3), 676-689.https://doi.org/10.1287/trsc.2019.0935
504
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
