77、经典VRP建模

import random as rd
import matplotlib.pyplot as plt
import math
from pyscipopt import Model, quicksum
from tool.Distance import calculate_distance

def points_generate(x_min, y_min, x_max, y_max, nums):

    points = [[0, (x_max+x_min)//2, (y_max+y_min)//2]]
    for i in range(1, nums):
        points.append([i, int((x_max-x_min) * rd.random() + x_min),
                     int((y_max-y_min) * rd.random() + y_min)])
    return points

def A_plot(points, A, x_min, y_min, x_max, y_max):
    scala = min(x_max - x_min, y_max - y_min)
    plt.xlim(x_min - 0.1 * scala, x_max + 0.1 * scala)
    plt.ylim(y_min - 0.1 * scala, y_max + 0.1 * scala)

    # 显示框大小
    plt.gcf().set_size_inches(8, 8)
    # 画点
    for i in range(len(points)):
        plt.text(points[i][1], points[i][2], str(points[i][0]), color='r', size=500 / len(A))

    # 画流量箭头
    arrow_params = {'length_includes_head': True, 'shape': 'right',
                    'head_starts_at_zero': True}
    # 箭头端点空出比例
    scala_edge = 0.1
    # 字符在边的比例
    scala_text = 0.3

    for i in range(len(points)):
        for j in range(len(points)):
            x, y, dx, dy = points[i][1], points[i][2], points[j][1] - points[i][1], points[j][2] - points[i][2]
            if (dx == 0 and dy == 0) or A[i][j] == 0:
                continue
            ddx, ddy = -0.01 * scala * dy / (math.sqrt(dx * dx + dy * dy)), 0.01 * scala * dx / (
                math.sqrt(dx * dx + dy * dy))
            ddx, ddy = 0, 0
            # plt.arrow(x + scala_edge * dx + ddx, y + scala_edge * dy + 0.01 * scala + ddy, (1 - 2 * scala_edge) * dx,
            #           (1 - 2 * scala_edge) * dy, width=0.01*scala, head_width=0.015*scala, **arrow_params)
            plt.arrow(x + scala_edge * dx + ddx, y + scala_edge * dy + 0.01 * scala + ddy, (1 - 2 * scala_edge) * dx,
                      (1 - 2 * scala_edge) * dy, width=0.01 * scala, head_width=0.0001, **arrow_params)
            plt.text(x + scala_text * dx + 4 * ddx, y + scala_text * dy + 4 * ddy, str(A[i][j]))

    plt.show()

#附：列出L的所有真子集
def subset(L):
    LL = []
    n = len(L)
    for i in range(1, pow(2, n)-1):
        L1 = []
        temp = i
        for j in range(n):
            if temp%2 == 1:
                L1.append(L[j])
            temp = temp//2
        LL.append(L1)
    return LL

def VRP_IP_directed(Dis, n, K, time_limit):

    model = Model("VRP_IP_directed")

    A = [[model.addVar(vtype="B", name="A[%s,%s]" % (i, j)) for j in range(n)] for i in range(n)]

    # 以总成本最小为目标
    model.setObjective(quicksum(Dis[i][j] * A[i][j] for i in range(n) for j in range(n)), "minimize")

    #几辆车的约束
    model.addCons(quicksum(A[0][j] for j in range(1, n)) == K)
    model.addCons(quicksum(A[j][0] for j in range(1, n)) == K)

    #其余一度约束
    for i in range(1, n):
        model.addCons(quicksum(A[i][j] for j in range(n) if i != j) == 1)
        model.addCons(quicksum(A[j][i] for j in range(n) if i != j) == 1)

    #无子环约束
    SS = subset([i for i in range(n)])
    for S in SS:
        model.addCons(quicksum(A[i][j] for i in S for j in range(n) if j not in S) >= 1)

    # 设置求解时间
    model.setRealParam("limits/time", time_limit)
    model.optimize()
    print("\ngap:", model.getGap())

    A1 = [[round(model.getVal(A[i][j])) for j in range(n)] for i in range(n)]

    return A1

def price(A, Dis):
    n = len(A)
    return sum([Dis[i][j] * A[i][j] for i in range(n) for j in range(n)])

if __name__ == "__main__":

    #传统VRP的定义，只限制车数，最小化运输距离
    #点的范围及个数
    x_min, y_min, x_max, y_max, nums = 0,0,100,100,15

    #点集生成
    points = points_generate(x_min, y_min, x_max, y_max, nums)
    Dis = [[calculate_distance(points[i][1], points[i][2], points[j][1], points[j][2], type) for j in range(nums)] for i in range(nums)]
    print(points)

    #车数
    K = 5

    A = VRP_IP_directed(Dis, len(Dis), K, 100)

    print("A:", price(A, Dis))

    A1 = [[1 if i < j and (A[i][j] == 1 or A[j][i] == 1) else 0 for j in range(len(points))] for i in range(len(points))]

    A_plot(points, A1, x_min, y_min, x_max, y_max)