【COPT】

COPT】"/>

【COPT】

Python调用COPT : 最短路问题及其对偶的探讨

• Python调用COPT：Shortest Path Problem及其对偶问题的一些探讨
• • Shortest Path Problem ： The Model
  • Python调用COPT求解SPP
  • SPP的对偶问题
  • Python调用COPT求解SPP的对偶问题
• SPP模型学术界的标准写法及其求解
• • SPP模型学术界的标准写法
  • SPP对偶问题的直观理解
• SPP模型学术界的标准形式的求解
• 参考文献

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

文中涉及到参考文献的标注方法：

Python调用COPT：Shortest Path Problem及其对偶问题的一些探讨

最短路问题 (Shortest Path Problem， SPP)是一类非常经典的问题。基本的SPP不是NP-hard,可以用Dijkstra等算法在多项式时间内求解到最优解。今天我们不探讨这些求解算法，仅探讨用求解器COPT和Python来求解这个问题。

我们首先来看一个例子网络：

s为起点，t为终点。

注意到边 b → a b \rightarrow a b→a的权重是-10，那么在这个问题中就存在负环，例如 b → a → c → b b \rightarrow a \rightarrow c \rightarrow b b→a→c→b，这个环的总权重为-5。这种情况下，对偶问题是无解的。用这个例子目的就是稍微拓展一下这一点。

Shortest Path Problem ： The Model

首先，我们将Shortest Path Problem建模成一个Integer Programming模型，如下：

min ⁡ ∑ ( i , j ) ∈ A x i j d i j ∑ ( s , j ) ∈ A x s j = 1 ∑ ( j , t ) ∈ A x j t = 1 ∑ ( i , j ) ∈ A x i j − ∑ ( j , k ) ∈ A x j k = 0 , ∀ j ∈ V \ { s , t } x i j ∈ { 0 , 1 } , ∀ ( i , j ) ∈ A \begin{aligned} \min &\sum_{\left( i,j \right) \in A}{x_{ij}d_{ij}} \\ &\sum_{\left( s,j \right) \in A}{x_{sj}}=1 \\ &\sum_{\left( j,t \right) \in A}{x_{jt}}=1 \\ &\sum_{\left( i,j \right) \in A}{x_{ij}}-\sum_{\left( j,k \right) \in A}{x_{jk}}=0, \quad \forall j\in V\backslash \left\{ s,t \right\} \\ &x_{ij}\in \left\{ 0,1 \right\} , \quad\forall \left( i,j \right) \in A \end{aligned} min​(i,j)∈A∑​xij​dij​(s,j)∈A∑​xsj​=1(j,t)∈A∑​xjt​=1(i,j)∈A∑​xij​−(j,k)∈A∑​xjk​=0,∀j∈V\{s,t}xij​∈{0,1},∀(i,j)∈A​

Python调用COPT求解SPP

下面我们用Python调用COPT，并用上图中的算例进行建模求解。

首先，导入模块，并将上述网络转换成字典格式数据

from coptpy import * Nodes = ['s', 'a', 'b', 'c', 't'] Arcs = {('s','a'): 5 ,('s','b'): 8,('a','c'): 2,('b','a'): -10,('c','b'): 3,('b','t'): 4,('c','t'): 3}
Arcs

然后建模

# Create COPT environment 
env = Envr()# create model 
model = env.createModel("SPP model") # add decision variables 
x = {} 
for key in Arcs.keys():x[key] = model.addVar(vtype=COPT.BINARY , name= 'x_' + key[0] + ',' + key[1] ) # add objective function
obj = LinExpr() 
for key in Arcs.keys():obj.addTerms(x[key], Arcs[key])
model.setObjective(obj, COPT.MINIMIZE)   # constraint1 1 and constraint 2  
lhs_1 = LinExpr()
lhs_2 = LinExpr()
for key in Arcs.keys():if(key[0] == 's'):lhs_1.addTerms(x[key], 1) elif(key[1] == 't'):lhs_2.addTerms(x[key], 1)
model.addConstr(lhs_1 == 1, name = 'start flow')
model.addConstr(lhs_2 == 1, name = 'end flow') # constraints 3
for node in Nodes:if(node != 's' and node != 't'):lhs = LinExpr()for key in Arcs.keys():if(key[1] == node):lhs.addTerms(x[key], 1)elif(key[0] == node):lhs.addTerms(x[key], -1)model.addConstr(lhs == 0, name = 'flow conservation')  model.write('model_spp.lp')
model.solve()print('objective value:', model.objVal) 
print('solution:') 
for var in model.getVars(): if(var.x > 0):print(var.getName(), ' = ', var.x)  

求解结果如下

objective value:Sat Feb 26 07:48:58 2022 [ Info] read client config: /opt/copt30/client.ini
Sat Feb 26 07:48:58 2022 [ Info] initialize floating client
Cardinal Optimizer v3.0.1. Build date Oct 12 2021
Copyright Cardinal Operations 2021. All Rights Reserved3.0
solution:
x_s,b  =  1.0
x_a,c  =  1.0
x_b,a  =  1.0
x_c,t  =  1.0

可以看到，最短路为 s → b → a → c → t s \rightarrow b \rightarrow a \rightarrow c \rightarrow t s→b→a→c→t，路径长度为3.

SPP的对偶问题

下面我们写出这个问题的对偶问题，我们对每个约束引入对偶变量 π i \pi_i πi​，如下

min ⁡ ∑ ( i , j ) ∈ A x i j d i j ∑ ( s , j ) ∈ A x s j = 1 → π s ∑ ( j , t ) ∈ A x j t = 1 → π t ∑ ( i , j ) ∈ A x i j − ∑ ( j , k ) ∈ A x j k = 0 , ∀ j ∈ V \ { s , t } → π j x i j ∈ { 0 , 1 } , ∀ ( i , j ) ∈ A \begin{aligned} \min &\sum_{\left( i,j \right) \in A}{x_{ij}d_{ij}} \\ &\sum_{\left( s,j \right) \in A}{x_{sj}}=1 \qquad \rightarrow \color{red} \quad \pi_s \\ &\sum_{\left( j,t \right) \in A}{x_{jt}}=1 \qquad \rightarrow \color{red} \quad \pi_t \\ &\sum_{\left( i,j \right) \in A}{x_{ij}}-\sum_{\left( j,k \right) \in A}{x_{jk}}=0, \quad \forall j\in V\backslash \left\{ s,t \right\} \qquad \rightarrow \color{red} \quad \pi_j \\ &x_{ij}\in \left\{ 0,1 \right\} , \quad\forall \left( i,j \right) \in A \end{aligned} min​(i,j)∈A∑​xij​dij​(s,j)∈A∑​xsj​=1→πs​(j,t)∈A∑​xjt​=1→πt​(i,j)∈A∑​xij​−(j,k)∈A∑​xjk​=0,∀j∈V\{s,t}→πj​xij​∈{0,1},∀(i,j)∈A​

然后SPP的对偶问题变成

max ⁡ π s + π t π s − π j ⩽ d s j , ∀ j ∈ { a , b } π j + π t ⩽ d j t , ∀ j ∈ { b , c } π i − π j ⩽ d i j , ∀ i , j ∈ A , i , j ∈ V \ { s , t } π free \begin{aligned} \max \quad &\pi _s+\pi _t \\ &\pi _s-\pi _j\leqslant d_{sj}, \qquad \forall j\in \left\{ a,b \right\} \\ &\pi _j+\pi _t\leqslant d_{jt}, \qquad \forall j\in \left\{ b,c \right\} \\ &\pi _i-\pi _j\leqslant d_{ij}, \qquad \forall i,j\in A, i,j \in V\backslash \left\{ s,t \right\} \\ &\pi \,\,\text{free} \end{aligned} max​πs​+πt​πs​−πj​⩽dsj​,∀j∈{a,b}πj​+πt​⩽djt​,∀j∈{b,c}πi​−πj​⩽dij​,∀i,j∈A,i,j∈V\{s,t}πfree​

Python调用COPT求解SPP的对偶问题

下面调用COPT进行求解

from coptpy import * env = Envr()model = env.createModel('dual') pi_a = model.addVar(lb=-1000, ub=1000, vtype=COPT.CONTINUOUS, name= "pi_a") 
pi_b = model.addVar(lb=-1000, ub=1000, vtype=COPT.CONTINUOUS, name= "pi_b") 
pi_c = model.addVar(lb=-1000, ub=1000, vtype=COPT.CONTINUOUS, name= "pi_c") 
pi_s = model.addVar(lb=-1000, ub=1000, vtype=COPT.CONTINUOUS, name= "pi_s")  
pi_t = model.addVar(lb=-1000, ub=1000, vtype=COPT.CONTINUOUS, name= "pi_t") obj = LinExpr() 
obj.addTerms(pi_s, 1)
obj.addTerms(pi_t, 1)
model.setObjective(obj, COPT.MAXIMIZE) # lhs relation , rhs 
model.addConstr(pi_s - pi_a <= 5)
model.addConstr(pi_s - pi_b <= 8)
model.addConstr(pi_t + pi_c <= 3)
model.addConstr(pi_t + pi_b <= 4) 
model.addConstr(pi_b - pi_a <= -10)    
model.addConstr(pi_a - pi_c <= 2)
model.addConstr(pi_c - pi_b <= 3)model.write('model_dual.mps')
model.solve()print('objective value:', model.objVal) 
print('solution:')  
for var in model.getVars():print(var.getName(), '\t', var.x) 

求解之后，结果为无解

Solving finished
Status: Infeasible  Objective: -  Iterations: 5  Time: 0.00s
Solution is not available

让我们修改边 b → a b \rightarrow a b→a的长度为10，也就是把下面的约束改成

# model.addConstr(pi_b - pi_a <= -10)
model.addConstr(pi_b - pi_a <= 10)

这样就有解了

objective value: 10.0
solution:
pi_a 	 -995.0
pi_b  -998.0
pi_c 	 -997.0
pi_s 	 -990.0
pi_t 	 1000.0

可以看到，最短路为 s → a → c → t s \rightarrow a \rightarrow c \rightarrow t s→a→c→t，路径长度为10.

SPP模型学术界的标准写法及其求解

SPP模型学术界的标准写法

上面的写法是比较适合初学者的，学者们一般不这么写，下面我们来整理一下比较专业的写法。主要参考文献¹ 和²。
我们将SPP的模型可以写为：
max ⁡ ∑ e ∈ A d e x e ∑ e ∈ o u t ( i ) x e − ∑ e ∈ i n ( i ) x e = b i , ∀ i ∈ V x e ∈ { 0 , 1 } , ∀ e ∈ A \begin{aligned} \max &\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}=b_i, \qquad \forall i\in V \\ &x_e\in \left\{ 0, 1 \right\} , \qquad \forall e\in A \end{aligned} max​e∈A∑​de​xe​e∈out(i)∑​xe​−e∈in(i)∑​xe​=bi​,∀i∈Vxe​∈{0,1},∀e∈A​
其中当 i = s i = s i=s时， b i = 1 b_i = 1 bi​=1, 当 i = t i=t i=t时， b i = − 1 b_i = -1 bi​=−1, 否则 b i = 0 b_i = 0 bi​=0。

当然，我这里写out, in什么的纯粹为了看着直观，其实也不是很直观，论文里有的是FS，BS。

也有写成下面形式的

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 x e ∈ { 0 , 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} \\ &x_e\in \left\{ 0,1 \right\} ,\qquad \forall e\in A \end{aligned} max​e∈A∑​de​xe​e∈out(i)∑​xe​−e∈in(i)∑​xe​=⎩⎪⎨⎪⎧​1,−1,0,​ifi=sifi=telse​xe​∈{0,1},∀e∈A​

看个人喜好了。

由于SPP具有整数最优解特性，也就是将决策变量 x e x_e xe​松弛成 0 ⩽ x e ⩽ 1 0 \leqslant x_e \leqslant 1 0⩽xe​⩽1,依然总是存在整数最优解.
因此，SPP等价于下面的LP。
min ⁡ ∑ e ∈ A d e x e ∑ e ∈ o u t ( i ) x e − ∑ e ∈ i n ( i ) x e = b i , ∀ i ∈ V 0 ⩽ x e ⩽ 1 , ∀ e ∈ A \begin{aligned} \min &\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}=b_i, \qquad \forall i\in V \\ &0 \leqslant x_e \leqslant 1 , \qquad \forall e\in A \end{aligned} min​e∈A∑​de​xe​e∈out(i)∑​xe​−e∈in(i)∑​xe​=bi​,∀i∈V0⩽xe​⩽1,∀e∈A​

相应的，对偶问题则变成
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πi​free​

到这一步其实已经结束了。我们可以做一个比较直观的理解(这个理解理论上有缺陷):

π i \pi_i πi​其实就是从 s s s点出发，到达点 j j j的收益的一个衡量。当然，数值上并不是这样的，之后我们做一步操作，就可以让数值上也相等了。

首先把所有对偶变量取个相反数，等价转化为
max ⁡ π t − π s π j − π i ⩽ d i j , ∀ ( i , j ) ∈ A π i free \begin{aligned} \max \qquad &\pi _t-\pi _s \\ &\pi _j-\pi _i \leqslant d_{ij}, \qquad \forall \left( i,j \right) \in A \\ &\pi _i\,\,\,\,\text{free} \end{aligned} max​πt​−πs​πj​−πi​⩽dij​,∀(i,j)∈Aπi​free​

因为 π i \pi_i πi​是free的，因此是等价的。

由于 π i \pi_i πi​是无约束的，而我们只是最大化 π s − π t \pi_s-\pi_t πs​−πt​，因此我们可以直接把 π s \pi_s πs​设置成0，也就是约束 π s = 0 \pi_s = 0 πs​=0，这样并不会改变最优解。OK，我们仍然保持目标函数为 max ⁡ \max max，那我们就可以把模型等价转化为

max ⁡ π t π j − π i ⩽ d i j , ∀ ( i , j ) ∈ A π i free , i ≠ s , π s = 0 \begin{aligned} \max \qquad &\pi _t \\ &\pi _j-\pi _i \leqslant d_{ij}, \qquad \forall \left( i,j \right) \in A \\ &\pi _i\,\,\,\,\text{free}, i \ne s, \pi _s = 0 \end{aligned} max​πt​πj​−πi​⩽dij​,∀(i,j)∈Aπi​free,i​=s,πs​=0​
这样我们就能保证所有的对偶变量 π i \pi_i πi​都是非负的，而且这个模型与原模型完全等价，而且问题复杂度相同。

SPP对偶问题的直观理解

依据这个模型，我们就可以对对问题进行一个比较深入的直观理解了。

就像生产计划问题的对偶问题，可以按照出租的逻辑来理解一样，SPP问题的原问题是一个客户要在网络上从起点到终点走一条最短路，花销最少。那SPP的对偶问题就可以理解为，我是网络的运营者，我需要在每个结点设置收费金额（或者惩罚金额），我设置的惩罚金额必须要使得我的收益在对方走最短路的情况下最大化。当然，对方不走最短路，我的收益会更大，这个惩罚，在一条边连接的两个节点上，差值不能超过这条边的权重。

这里，我们戏称运营者为地主。地主为了让平民走得最短，就得拼命的加惩罚，狠狠收费。怎么衡量地主收费收的是不是狠呢？下面给出一些直观理解的分享(并不完全严谨，仅供参考，帮助理解)。

一条边 ( i , j ) (i,j) (i,j)的收费上限是 d i j d_{ij} dij​。地主收费不狠，那就是 π j − π i < < d i j \pi_j - \pi_i << d_{ij} πj​−πi​<<dij​，说明他还有点良心。余量越多，他的收益就越小，对应的，原问题相应的解就越远离最短路。收的越狠，接近 π j − π i ≈ d i j \pi_j - \pi_i \approx d_{ij} πj​−πi​≈dij​了，原问题（平民）就不敢瞎浪了，乖乖走到最短路，毕竟财力限制了他的想象力。当 π j − π i = d i j \pi_j - \pi_i =d_{ij} πj​−πi​=dij​，这就使得原问题走到最短路。

改进的对偶问题的解，也就是原问题的对偶变量，正好也符合上面的理解。

SPP模型学术界的标准形式的求解

按照修改的模型，我们再将今天的例子求解一下，约束按照修改的模型进行构建。

from coptpy import * env = Envr()model = env.createModel('dual') pi_a = model.addVar(lb=0, ub=1000, vtype=COPT.CONTINUOUS, name= "pi_a") 
pi_b = model.addVar(lb=0, ub=1000, vtype=COPT.CONTINUOUS, name= "pi_b") 
pi_c = model.addVar(lb=0, ub=1000, vtype=COPT.CONTINUOUS, name= "pi_c") 
pi_s = model.addVar(lb=0, ub=1000, vtype=COPT.CONTINUOUS, name= "pi_s")  
pi_t = model.addVar(lb=0, ub=1000, vtype=COPT.CONTINUOUS, name= "pi_t") obj = LinExpr() 
obj.addTerms(pi_s, 1)
obj.addTerms(pi_t, 1)
model.setObjective(obj, COPT.MAXIMIZE) # lhs relation , rhs 
model.addConstr(pi_s - pi_a <= 5)
model.addConstr(pi_s - pi_b <= 8)
model.addConstr(pi_t + pi_c <= 3)
model.addConstr(pi_t + pi_b <= 4) 
model.addConstr(pi_b - pi_a <= 10)    # 如果是-10，就会不可行
model.addConstr(pi_a - pi_c <= 2)
model.addConstr(pi_c - pi_b <= 3)model.write('model_dual.mps')
model.solve()print('objective value:', model.objVal) 
print('solution:')  
for var in model.getVars():print(var.getName(), '\t', var.x) 

求解结果如下

Solved in 2 iterations and 0.01 seconds
Optimal objective  1.000000000e+01
10.0
pi_a 	 5.0
pi_b 	 6.0
pi_c 	 7.0
pi_s 	 0.0
pi_t 	 10.0

上述结果跟前面的理解是相同的。 π i , ∀ i ∈ V \pi_i, \forall i \in V πi​,∀i∈V可以看做是在点 i i i处地主能够收到的过路费的最大值。在终点能够收到的过路费的最大值 π t \pi_t πt​就是最终的目标函数。而过路费必须满足决策者提出的 π j − π i ⩽ d i j \pi_j - \pi_i \leqslant d_{ij} πj​−πi​⩽dij​的反垄断要求。

然后平民就乖乖走出了最短路。这问题就这么愉快的解决了。

参考文献

- [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-652. - [2]: Cappanera, P., & Scaparra, M. P. (2011). Optimal allocation of protective resources in shortest-path networks. Transportation Science, 45(1), 64-80.

---

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

公众号往期推文如下