python minimize_Python数学规划案例三：最短路径

Python数学规划案例三：最短路径

在图论中，最短路径问题(the shortest path problem)是在图中的两个顶点(或节点)之间寻找路径，使得构成路径的所有边的权重之和最小，权重通常为距离成本。

数学模型

数据集

来源：

http://people.brunel.ac.uk/~mastjjb/jeb/orlib/rcspinfo.html

具体说明见问题二数据集描述，所需数据为节点(N)以及弧成本(C_ij)。

Python代码

0) 载入packages

from docplex.mp.model import Model

import numpy as np

1) 把算例数据集文件下载下来，本文的例子使用rcsp5.txt。

2) 生成class SPP，指定读写文件目录。

class SPP:

def __init__(self):

self.mDir = 'D:\\PY_CPLEX\\SPP\\'

self.mDir_input = 'i\\'

self.mDir_output = 'o\\'

self.mData = {}

# mData:{'vertices_num','arcs_num','resources_num','resources_limit','vertices_consume','arcs_cost_consume',}

def read_data(self):

......

def solve_by_cplex(self):

......

def experiment(self):

......

if __name__ == '__main__':

......

3) 定义读数据的函数read_data。

def read_data(self):

ls_fn = self.mDir + self.mDir_input

filename = 'rcsp5.txt'

with open(ls_fn + filename, 'r') as f: # 读取数据

data = f.readlines()

data_list = []

for data_str in data:

data_list += data_str.split() # data_list = ['100','990','10','0','0',...]

k = 0

vertices_num = int(data_list[k]) # (1)节点数量 100

k += 1

arcs_num = int(data_list[k]) # (2)弧的数量是 990

k += 1

resources_num = int(data_list[k]) # (3)资源数量是 10

k += 1

Lk = []

Uk = []

for i in range(2): #i=[0,1]

for j in range(resources_num):

if i == 0:

Lk.append(int(data_list[k])) #Lk = ['0','0','0',...]资源下限

k += 1

else:

Uk.append(int(data_list[k]))

k += 1 # Uk = ['178','170','167,'121',...]资源上限

Qik = [[0 for j in range(resources_num)] for i in range(vertices_num)]

for i in range(vertices_num):

for j in range(resources_num):

Qik[i][j] = int(data_list[k])

k += 1 # 经过节点i消耗的资源k的数量，未经过的节点资源消耗数量设置为0

Rijk = [[[99999 for o in range(resources_num)] for j in range(vertices_num)] for i in range(vertices_num)]

# 经过弧(i,j)消耗资源k的数量，不通的弧设置为99999

cij = [[99999 for j in range(vertices_num)] for i in range(vertices_num)]

# 经过弧(i,j)需要的成本，不通的弧设置为99999

for i in range(arcs_num):

p1 = int(data_list[k]) - 1

k += 1

p2 = int(data_list[k]) - 1

k += 1

cij[p1][p2] = int(data_list[k])

k += 1

for j in range(resources_num):

Rijk[p1][p2][j] = int(data_list[k]) #Rijk[0][16][0]=13, Rijk[0][16][1]=14

k += 1

self.mData['vertices_num'] = vertices_num # 节点数量

self.mData['arcs_num'] = arcs_num # 弧的数量

self.mData['cij'] = cij # 经过弧ij所需要的成本

4) 定义求解数学规划模型的函数solvebycplex。注意其中涉及到的cplex的函数或属性，例如binaryvarmatrix，minimize，sum，add_constraint，add_constraints，solve，solution_value。

def solve_by_cplex(self):

# load data from properties

# STEP 1 Set

N = self.mData['vertices_num'] #节点的集合

# STEP 2 Parameter

cij = self.mData['cij']

# STEP 3 model

m = Model("spp")

# STEP 4 Variable

x = m.binary_var_matrix(N, N, name="x")

# xij = m.addVars(N, N, vtype=GRB.BINARY, name="xij")

# STEP 5 Objective

m.minimize(m.sum(cij[i][j] * x[i, j] for i in range(N) for j in range(N)))

# STEP 6 Constraint

#m.add_constraints(m.sum(x[i, j] for i in range(N)) == 1 for j in range(N))

m.add_constraint(m.sum(x[0, j] for j in range(N)) == 1) # (1)起点的出度为1,python索引从0开始

m.add_constraint(m.sum(x[i, N-1] for i in range(N)) == 1) # (2)终点的入度为1,最后一个点的索引为N-1

m.add_constraints(m.sum(x[i, j] for i in range(N)) == m.sum(x[j, i] for i in range(N)) for j in range(1,N-1))

# (3)除起点和终点外，每个节点的出度等于入度

# STEP 7 solve

s = m.solve(log_output=True)

# STEP 8 save result to file

filename = self.mDir + self.mDir_output + 'spp5_result.txt'

file_handle = open(filename, 'w')

file_handle.write('The minimal cost is %g\n' % m.objective_value)

file_handle.write('Routing:\n[From, To]\n')

mat = np.zeros([N, N])

for i in range(N):

for j in range(N):

if x[i, j].solution_value > 0.9:

mat[i][j] = x[i, j]

[a, b] = np.where(mat == 1)

sqe1 = [0]

for i in range(1, len(a)+1):

id = np.where(a == sqe1[i-1])

sqe1.append(b[id])

file_handle.write("[%d, %d]\n" % (sqe1[i-1], sqe1[i]))

print("From_vertex_%d, To_vertex_%d]\n" % (sqe1[i-1], sqe1[i]))

file_handle.close()

5) 定义函数experiment执行读数据和求解。

def experiment(self):

self.read_data()

self.solve_by_cplex()

6) 编写main部分代码。

if __name__ == '__main__':

a = SPP()

a.experiment()

完整的代码

from docplex.mp.model import Model

import numpy as np

class SPP:

def __init__(self):

self.mDir = 'D:\\PY_CPLEX\\SPP\\'

self.mDir_input = 'i\\'

self.mDir_output = 'o\\'

self.mData = {}

# mData:{'vertices_num','arcs_num','resources_num','resources_limit','vertices_consume','arcs_cost_consume',}

def read_data(self):

ls_fn = self.mDir + self.mDir_input

filename = 'rcsp5.txt'

with open(ls_fn + filename, 'r') as f: # 读取数据

data = f.readlines()

data_list = []

for data_str in data:

data_list += data_str.split() # data_list = ['100','990','10','0','0',...]

k = 0

vertices_num = int(data_list[k]) # (1)节点数量 100

k += 1

arcs_num = int(data_list[k]) # (2)弧的数量是 990

k += 1

resources_num = int(data_list[k]) # (3)资源数量是 10

k += 1

Lk = []

Uk = []

for i in range(2): #i=[0,1]

for j in range(resources_num):

if i == 0:

Lk.append(int(data_list[k])) #Lk = ['0','0','0',...]资源下限

k += 1

else:

Uk.append(int(data_list[k]))

k += 1 # Uk = ['178','170','167,'121',...]资源上限

Qik = [[0 for j in range(resources_num)] for i in range(vertices_num)]

for i in range(vertices_num):

for j in range(resources_num):

Qik[i][j] = int(data_list[k])

k += 1 # 经过节点i消耗的资源k的数量，未经过的节点资源消耗数量设置为0

Rijk = [[[99999 for o in range(resources_num)] for j in range(vertices_num)] for i in range(vertices_num)]

# 经过弧(i,j)消耗资源k的数量，不通的弧设置为99999

cij = [[99999 for j in range(vertices_num)] for i in range(vertices_num)]

# 经过弧(i,j)需要的成本，不通的弧设置为99999

for i in range(arcs_num):

p1 = int(data_list[k]) - 1

k += 1

p2 = int(data_list[k]) - 1

k += 1

cij[p1][p2] = int(data_list[k])

k += 1

for j in range(resources_num):

Rijk[p1][p2][j] = int(data_list[k]) #Rijk[0][16][0]=13, Rijk[0][16][1]=14

k += 1

self.mData['vertices_num'] = vertices_num # 节点数量

self.mData['arcs_num'] = arcs_num # 弧的数量

self.mData['cij'] = cij # 经过弧ij所需要的成本

def solve_by_cplex(self):

# load data from properties

# STEP 1 Set

N = self.mData['vertices_num'] #节点的集合

# STEP 2 Parameter

cij = self.mData['cij']

# STEP 3 model

m = Model("spp")

# STEP 4 Variable

x = m.binary_var_matrix(N, N, name="x")

# xij = m.addVars(N, N, vtype=GRB.BINARY, name="xij")

# STEP 5 Objective

m.minimize(m.sum(cij[i][j] * x[i, j] for i in range(N) for j in range(N)))

# STEP 6 Constraint

#m.add_constraints(m.sum(x[i, j] for i in range(N)) == 1 for j in range(N))

m.add_constraint(m.sum(x[0, j] for j in range(N)) == 1) # (1)起点的出度为1,python索引从0开始

m.add_constraint(m.sum(x[i, N-1] for i in range(N)) == 1) # (2)终点的入度为1,最后一个点的索引为N-1

m.add_constraints(m.sum(x[i, j] for i in range(N)) == m.sum(x[j, i] for i in range(N)) for j in range(1,N-1))

# (3)除起点和终点外，每个节点的出度等于入度

# STEP 7 solve

s = m.solve(log_output=True)

# STEP 8 save result to file

filename = self.mDir + self.mDir_output + 'spp5_result.txt'

file_handle = open(filename, 'w')

file_handle.write('The minimal cost is %g\n' % m.objective_value)

file_handle.write('Routing:\n[From, To]\n')

mat = np.zeros([N, N])

for i in range(N):

for j in range(N):

if x[i, j].solution_value > 0.9:

mat[i][j] = x[i, j]

[a, b] = np.where(mat == 1)

sqe1 = [0]

for i in range(1, len(a)+1):

id = np.where(a == sqe1[i-1])

sqe1.append(b[id])

file_handle.write("[%d, %d]\n" % (sqe1[i-1], sqe1[i]))

print("From_vertex_%d, To_vertex_%d]\n" % (sqe1[i-1], sqe1[i]))

file_handle.close()

def experiment(self):

self.read_data()

self.solve_by_cplex()

if __name__ == '__main__':

a = SPP()

a.experiment()