运筹学修炼日记：如何优雅地写出大规模线性规划的对偶问题

运筹学修炼日记：如何优雅地写出大规模线性规划的对偶
- 最短路问题
- 多商品流问题`Multicommodity Network Flow Problem`
借助`Excel`和`具体小算例`写出大规模LP的对偶
- Dual Problem :Shortest Path Problem(最短路问题)
- - 小算例
- `Excel`+`小算例`写出`SPP`的对偶问题
- - 将`Excel`中的`Dual tabular`转化成公式形式
- Dual Problem :Multicommodity Network Flow Problem(最短路问题)
- - 具体算例
- `Excel`+`小算例`写出`MNF`的原模型
- - 将`Excel`中的`Dual tabular`转化成公式形式
- Python调用Gurobi求解Multicommodity Network Flow Problem （仅原问题）
- 后记

欢迎关注我们的微信公众号 运小筹

运筹学修炼日记：如何优雅地写出大规模线性规划的对偶

对偶理论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 mathArcs = {'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]}
ArcsNodes = [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

欢迎关注我们的微信公众号 运小筹

公众号往期推文如下