CPLEX Python【1】--DoCplex建模解CVRP问题

学习笔记而已，看的b站搬运
https://www.bilibili.com/video/BV1ab411G79M
CVRP指载重量有限的车辆运输问题，从起点出发，到n个点获取物资，然后运回起点，要求运输路径总和最短
主要把视频里面每一句代码的意义搞清楚(注释)，方便自己复现学习与对照使用

import cplex
import docplex
import numpy as np
rnd = np.random
rnd.seed(0)


n = 10                                #10个装载点
Q = 15                                #15个装载容量
N = [i for i in range(1, n+1)]        #10个运输点的各自编号
V = [0] + N                           #所有顶点编号，以0为出发顶点
q = {i:rnd.randint(1, 10) for i in N} #每个装载点拥有的货物量

loc_x = rnd.rand(len(V))*200
loc_y = rnd.rand(len(V))*100       #为每一个顶点在平面坐标系上随机一个位置，之后点间距将以这个为基础计算


import matplotlib.pyplot as plt
plt.scatter(loc_x[1:],loc_y[1:],c='b')
for i in N:
    plt.annotate("$q_%d=%d$"%(i, q[i]),(loc_x[i] + 2, loc_y[i]))
plt.scatter(loc_x[0],loc_y[0],c='r',marker='s')#显示随机出来的各个点，其中红色的定义为起点

生成点如图：
这里点数是15，但其实作为np问题15个点要算很久很久（算了100多秒才得到最优解），所以代码里换成了10

设置点间距
视频中对点的批量处理基本上都用一些简单的语法糖
我为了方便学习，都使用for循环实现

A = []
c = {}
for i in V:
    for j in V:
        if (i!=j):
            A.append((i, j))  #创建所有点对
for i,j in A:
    c[i,j] = np.hypot(loc_x[i] - loc_x[j], loc_y[i] - loc_y[j]) #定义每个点对之间的距离为坐标系平面距离


from docplex.mp.model import Model


mdl = Model(name="Test")
x = mdl.binary_var_dict(A, name='x') #创建自变量集x，x(i,j)为 1 表示选择该条边
u = mdl.continuous_var_dict(N, ub=Q, name='u') #创建自变量集u，u(i)表示到i点后的车辆载重

添加约束限制，进行优化求解

mdl.minimize(mdl.sum(c[i,j]*x[i,j] for i,j in A))  #设定目标，让所有选择的边代价总和最小

#添加限制1，这一行限制是指从i出发的边只能选择一个
for i in N:
    mdl.add_constraint(mdl.sum(x[i,j] for j in V if j!=i) == 1)  #使用sum函数还是适合用用语法糖
#添加限制2，这一行限制是指到达j的边中只能选择一个
for j in N:
    mdl.add_constraint(mdl.sum(x[i,j] for i in V if j!=i) == 1)
#限制1、2的内循环区间为V，而外循环为N，是因为对于出发点不需要规定只能有一条边出去或进来

#indicator，条件限制。
#这里设置多条条件限制，每一条条件限制意义为
#当选择了x[i,j]这一条边时，就限制到j点的总载重量必须等于到i的加上j所拥有的物资量
#这里指定一下i、j都不能为0，因为u[0]\q[0] 没有定义 并且 没有现实意义
for i,j in A:
    if i!=0 and j!=0:
        mdl.add_indicator_constraint(mdl.indicator_constraint(x[i,j], u[i] + q[j] == u[j]))

#限制一下到达i点后的总载重量，显然任何情况下都应该大于i点本身拥有的物资量（只经过一次）
#更好的理解应该是给u设置初始值
mdl.add_constraints(u[i]>=q[i] for i in N)

solution = mdl.solve(log_output=True)

显示答案

used_edge = [(i,j) for i,j in A if x[(i,j)].solution_value > 0.9] #得到选用的各个边

import matplotlib.pyplot as plt
plt.scatter(loc_x[1:],loc_y[1:],c='b')
for i in N:
    plt.annotate("$q_%d=%d$"%(i, q[i]),(loc_x[i] + 2, loc_y[i]))
for i,j in used_edge:
    plt.plot([loc_x[i], loc_x[j]],[loc_y[i], loc_y[j]],c='k')
plt.scatter(loc_x[0],loc_y[0],c='r',marker='s')#显示随机出来的各个点，其中红色的定义为起点