37、带约束的选址路径问题

import random as rd
import matplotlib.pyplot as plt
import functools
import math
import copy
import time
from pyscipopt import Model, quicksum

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

    points = []
    for i in range(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 OD_generate(num_min, num_max, nums_w, nums_d, nums_s):
    ODs = [[0 for j in range(nums_w+nums_d+nums_s)] for i in range(nums_w+nums_d+nums_s)]

    for i in range(nums_w):
        for j in range(nums_w+nums_d, nums_w+nums_d+nums_s):
            ODs[i][j] = int((num_max-num_min) * rd.random() + num_min)

    return ODs

def constraint_generate(nums_w, nums_d, nums_s, capacity_low, capacity_up):

    constraint_w_to_d = [[2 for j in range(nums_d)] for i in range(nums_w)]

    #假设0.1的概率有强制（1约束）约束，其余点0.5概率可能，0.5概率不可能
    # for i in range(nums_w):
    #     flag = rd.choice([j for j in range(10)])
    #     if flag == 0:
    #         idex_delv = rd.choice([j for j in range(nums_d)])
    #         for j in range(nums_d):
    #             constraint_w_to_d[i][j] = 0
    #         constraint_w_to_d[i][idex_delv] = 1
    #     else:
    #         for j in range(nums_d):
    #             constraint_w_to_d[i][j] = rd.choice([0,2])
    #         if sum(constraint_w_to_d[i]) == 0:
    #             idex_delv = rd.choice([j for j in range(nums_d)])
    #             constraint_w_to_d[i][idex_delv] = 2

    constraint_d_to_s = [[2 for j in range(nums_s)] for i in range(nums_d)]

    #假设0.1的概率有强制（1约束）约束，其余点0.5概率可能，0.5概率不可能
    # for j in range(nums_s):
    #     flag = rd.choice([i for i in range(10)])
    #     if flag == 0:
    #         idex_delv = rd.choice([i for i in range(nums_d)])
    #         for i in range(nums_d):
    #             constraint_d_to_s[i][j] = 0
    #         constraint_d_to_s[idex_delv][j] = 1
    #     else:
    #         for i in range(nums_d):
    #             constraint_d_to_s[i][j] = rd.choice([0,2])
    #         if sum([constraint_d_to_s[i][j] for i in range(nums_d)]) == 0:
    #             idex_delv = rd.choice([i for i in range(nums_d)])
    #             constraint_d_to_s[idex_delv][j] = 2

    capacity = [int((capacity_up - capacity_low) * rd.random() + capacity_up) for i in range(nums_d)]

    return constraint_w_to_d, constraint_d_to_s, capacity

def calculate_euclidean_distance(x1, y1, x2, y2):
    return math.sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))

def R_generate_mix_cost_constraint(points, nums_w, nums_d, nums_s, ODs, vehicles, constraint_w_to_d, constraint_d_to_s, capacity, time_limit):

    model = Model("Level_mix_cost_constraint")

    M = sum([sum(ODs[i]) for i in range(len(ODs))])*10
    print(M)
    R = [[[[model.addVar(vtype="B", name="R[%s,%s,%s,%s]" % (i, j, k, l)) for l in range(nums_s)] for k in range(nums_d)] for j in range(nums_d)] for i in range(nums_w)]

    bills_route_w_to_d = [[model.addVar(vtype="I", name="bills_route_w_to_d[%s,%s]" % (i, j)) for j in range(nums_d)] for i in range(nums_w)]
    bills_route_d_to_d = [[model.addVar(vtype="I", name="bills_route_d_to_d[%s,%s]" % (i, j)) for j in range(nums_d)] for i in range(nums_d)]
    bills_route_d_to_s = [[model.addVar(vtype="I", name="bills_route_d_to_s[%s,%s]" % (i, j)) for j in range(nums_s)] for i in range(nums_d)]
    nums_vehicles_w_to_d = [[[model.addVar(vtype="I", name="nums_vehicles_w_to_d[%s,%s,%s]" % (i, j, k)) for k in range(len(vehicles))] for j in range(nums_d)] for i in range(nums_w)]
    nums_vehicles_d_to_d = [[[model.addVar(vtype="I", name="nums_vehicles_d_to_d[%s,%s,%s]" % (i, j, k)) for k in range(len(vehicles))] for j in range(nums_d)] for i in range(nums_d)]
    nums_vehicles_d_to_s = [[[model.addVar(vtype="I", name="nums_vehicles_d_to_s[%s,%s,%s]" % (i, j, k)) for k in range(len(vehicles))] for j in range(nums_s)] for i in range(nums_d)]

    # 欧式距离和球面距离
    Dis = [[calculate_euclidean_distance(points[i][1], points[i][2], points[j][1], points[j][2]) for j in range(nums_w + nums_d + nums_s)] for i in range(nums_w + nums_d + nums_s)]

    # 以总成本最小为目标
    model.setObjective(quicksum((2 * Dis[i][j+nums_w] * vehicles[k][3] + vehicles[k][2]) * nums_vehicles_w_to_d[i][j][k] for i in range(nums_w) for j in range(nums_d) for k in range(len(vehicles))) + quicksum((2 * Dis[i+nums_w][j+nums_w] * vehicles[k][3] + vehicles[k][2]) * nums_vehicles_d_to_d[i][j][k] for i in range(nums_d) for j in range(nums_d) for k in range(len(vehicles)) if i != j) + quicksum((2 * Dis[i+nums_w][j+nums_w+nums_d] * vehicles[k][3] + vehicles[k][2]) * nums_vehicles_d_to_s[i][j][k] for i in range(nums_d) for j in range(nums_s) for k in range(len(vehicles))), "minimize")

    # 变量约束
    for i in range(nums_w):
        for j in range(nums_d):
            model.addCons(bills_route_w_to_d[i][j] >= 0)
            model.addCons(bills_route_w_to_d[i][j] <= M)
            for k in range(len(vehicles)):
                model.addCons(nums_vehicles_w_to_d[i][j][k] >= 0)
                model.addCons(nums_vehicles_w_to_d[i][j][k] <= M)

    for i in range(nums_d):
        for j in range(nums_d):
            model.addCons(bills_route_d_to_d[i][j] >= 0)
            model.addCons(bills_route_d_to_d[i][j] <= M)
            for k in range(len(vehicles)):
                model.addCons(nums_vehicles_d_to_d[i][j][k] >= 0)
                model.addCons(nums_vehicles_d_to_d[i][j][k] <= M)

    for i in range(nums_d):
        for j in range(nums_s):
            model.addCons(bills_route_d_to_s[i][j] >= 0)
            model.addCons(bills_route_d_to_s[i][j] <= M)
            for k in range(len(vehicles)):
                model.addCons(nums_vehicles_d_to_s[i][j][k] >= 0)
                model.addCons(nums_vehicles_d_to_s[i][j][k] <= M)

    for i in range(nums_w):
        for j in range(nums_d):
            for k in range(nums_d):
                for l in range(nums_s):
                    model.addCons(R[i][j][k][l] >= 0)
                    model.addCons(R[i][j][k][l] <= 1)

    #单条路线约束
    for i in range(nums_w):
        for l in range(nums_s):
            model.addCons(quicksum(R[i][j][k][l] for j in range(nums_d) for k in range(nums_d)) == 1)

    #路线约束
    for i in range(nums_w):
        flag = 0
        idex_delv = 0
        for j in range(nums_d):
            if constraint_w_to_d[i][j] == 1:
                flag = 1
                idex_delv = j
                break
        if flag == 1:
            for j in range(nums_d):
                if j != idex_delv:
                    model.addCons(bills_route_w_to_d[i][j] == 0)
        else:
            for j in range(nums_d):
                if constraint_w_to_d[i][j] == 0:
                    model.addCons(bills_route_w_to_d[i][j] == 0)

    for j in range(nums_s):
        flag = 0
        idex_delv = 0
        for i in range(nums_d):
            if constraint_d_to_s[i][j] == 1:
                flag = 1
                idex_delv = i
                break
        if flag == 1:
            for i in range(nums_d):
                if i != idex_delv:
                    model.addCons(bills_route_d_to_s[i][j] == 0)
        else:
            for i in range(nums_d):
                if constraint_d_to_s[i][j] == 0:
                    model.addCons(bills_route_d_to_s[i][j] == 0)

    #分拣能力约束
    for i in range(nums_d):
        model.addCons(quicksum(bills_route_d_to_d[i][j] for j in range(nums_d)) + quicksum(bills_route_d_to_s[i][j] for j in range(nums_s)) <= capacity[i])

    #运输单量平衡
    for i in range(nums_w):
        for j in range(nums_d):
            model.addCons(bills_route_w_to_d[i][j]-quicksum(ODs[i][l+nums_w+nums_d]*R[i][j][k][l] for k in range(nums_d) for l in range(nums_s)) == 0)

    for j in range(nums_d):
        for k in range(nums_d):
            model.addCons(bills_route_d_to_d[j][k]-quicksum(ODs[i][l+nums_w+nums_d]*R[i][j][k][l] for i in range(nums_w) for l in range(nums_s)) == 0)

    for k in range(nums_d):
        for l in range(nums_s):
            model.addCons(bills_route_d_to_s[k][l]-quicksum(ODs[i][l+nums_w+nums_d]*R[i][j][k][l] for i in range(nums_w) for j in range(nums_d)) == 0)

    #车辆运载约束
    for i in range(nums_w):
        for j in range(nums_d):
            model.addCons(bills_route_w_to_d[i][j] - quicksum(nums_vehicles_w_to_d[i][j][k]*vehicles[k][1] for k in range(len(vehicles))) <= 0)

    for i in range(nums_d):
        for j in range(nums_d):
            if i != j:
                model.addCons(bills_route_d_to_d[i][j] - quicksum(nums_vehicles_d_to_d[i][j][k]*vehicles[k][1] for k in range(len(vehicles))) <= 0)

    for i in range(nums_d):
        for j in range(nums_s):
            model.addCons(bills_route_d_to_s[i][j] - quicksum(nums_vehicles_d_to_s[i][j][k]*vehicles[k][1] for k in range(len(vehicles))) <= 0)

    # 设置求解时间
    model.setRealParam("limits/time", time_limit)

    model.optimize()
    print("\ngap:",model.getGap())

    # 拿结果
    R1 = [[[[round(model.getVal(R[i][j][k][l])) for l in range(nums_s)] for k in range(nums_d)] for j in range(nums_d)] for i in range(nums_w)]
    bills_route_w_to_d1 = [[round(model.getVal(bills_route_w_to_d[i][j])) for j in range(nums_d)] for i in range(nums_w)]
    bills_route_d_to_d1 = [[round(model.getVal(bills_route_d_to_d[i][j])) for j in range(nums_d)] for i in range(nums_d)]
    bills_route_d_to_s1 = [[round(model.getVal(bills_route_d_to_s[i][j])) for j in range(nums_s)] for i in range(nums_d)]

    ODs1 = [[0 for j in range(nums_w+nums_d+nums_s)] for i in range(nums_w+nums_d+nums_s)]

    #车辆运载约束
    for i in range(nums_w):
        for j in range(nums_d):
            ODs1[i][j+nums_w] = bills_route_w_to_d1[i][j]

    for i in range(nums_d):
        for j in range(nums_d):
            ODs1[i+nums_w][j+nums_w] = bills_route_d_to_d1[i][j]

    for i in range(nums_d):
        for j in range(nums_s):
            ODs1[i+nums_w][j+nums_w+nums_d] = bills_route_d_to_s1[i][j]

    return R1, ODs1

def OD_plot(points, ODs, 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=scala/10)

    #画流量箭头
    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 ODs[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))
            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.text(x + scala_text*dx + 4*ddx, y + scala_text*dy + 4*ddy, str(ODs[i][j]))

    plt.show()

def check_constraint(nums_w, nums_d, nums_s, bills_route_w_to_d, nums_vehicles_w_to_d, bills_route_d_to_d, nums_vehicles_d_to_d, bills_route_d_to_s, nums_vehicles_d_to_s,
                     constraint_w_to_d, constraint_d_to_s, capacity, R):

    print("check bills_route_w_to_d、nums_vehicles_w_to_d")

    for i in range(nums_w):
        for j in range(nums_d):
            if bills_route_w_to_d[i][j] < 0:
                return False
            for k in range(len(vehicles)):
                if nums_vehicles_w_to_d[i][j][k] < 0:
                    return False

    print("check bills_route_d_to_d、nums_vehicles_d_to_d")

    for i in range(nums_d):
        for j in range(nums_d):
            if bills_route_d_to_d[i][j] < 0:
                return False
            for k in range(len(vehicles)):
                if nums_vehicles_d_to_d[i][j][k] < 0:
                    return False

    print("check bills_route_d_to_s、nums_vehicles_d_to_s")

    for i in range(nums_d):
        for j in range(nums_s):
            if bills_route_d_to_s[i][j] < 0:
                return False
            for k in range(len(vehicles)):
                if nums_vehicles_d_to_s[i][j][k] < 0:
                    return False

    #单条路线约束
    print("check 单条路线约束")
    for i in range(nums_w):
        for l in range(nums_s):
            if sum([R[i][j][k][l] for j in range(nums_d) for k in range(nums_d)]) != 1:
                return False

    #路线约束
    print("check 路线约束")
    for i in range(nums_w):
        flag = 0
        idex_delv = 0
        for j in range(nums_d):
            if constraint_w_to_d[i][j] == 1:
                flag = 1
                idex_delv = j
                break
        if flag == 1:
            for j in range(nums_d):
                if j != idex_delv:
                    if bills_route_w_to_d[i][j] != 0:
                        return False
        else:
            for j in range(nums_d):
                if constraint_w_to_d[i][j] == 0:
                    if bills_route_w_to_d[i][j] != 0:
                        return False

    print("check 路线约束2")
    for j in range(nums_s):
        flag = 0
        idex_delv = 0
        for i in range(nums_d):
            if constraint_d_to_s[i][j] == 1:
                flag = 1
                idex_delv = i
                break
        if flag == 1:
            for i in range(nums_d):
                if i != idex_delv:
                    if bills_route_d_to_s[i][j] != 0:
                        return False
        else:
            for i in range(nums_d):
                if constraint_d_to_s[i][j] == 0:
                    if bills_route_d_to_s[i][j] != 0:
                        return False

    #分拣能力约束
    print("check 分拣能力约束")
    for i in range(nums_d):
        if sum([bills_route_d_to_d[i][j] for j in range(nums_d)]) + sum([bills_route_d_to_s[i][j] for j in range(nums_s)]) > capacity[i]:
            return False

    #运输单量平衡
    print("check 运输单量平衡")
    for i in range(nums_w):
        for j in range(nums_d):
            if bills_route_w_to_d[i][j]-sum([ODs[i][l+nums_w+nums_d]*R[i][j][k][l] for k in range(nums_d) for l in range(nums_s)]) != 0:
                return False

    for j in range(nums_d):
        for k in range(nums_d):
            if bills_route_d_to_d[j][k]-sum([ODs[i][l+nums_w+nums_d]*R[i][j][k][l] for i in range(nums_w) for l in range(nums_s)]) != 0:
                return False

    for k in range(nums_d):
        for l in range(nums_s):
            if bills_route_d_to_s[k][l]-sum([ODs[i][l+nums_w+nums_d]*R[i][j][k][l] for i in range(nums_w) for j in range(nums_d)]) != 0:
                return False

    #车辆运载约束
    print("check 车辆运载约束")
    for i in range(nums_w):
        for j in range(nums_d):
            if bills_route_w_to_d[i][j] - sum([nums_vehicles_w_to_d[i][j][k]*vehicles[k][1] for k in range(len(vehicles))]) > 0:
                return False

    for i in range(nums_d):
        for j in range(nums_d):
            if i != j:
                if bills_route_d_to_d[i][j] - sum([nums_vehicles_d_to_d[i][j][k]*vehicles[k][1] for k in range(len(vehicles))]) > 0:
                    return False

    for i in range(nums_d):
        for j in range(nums_s):
            if bills_route_d_to_s[i][j] - sum([nums_vehicles_d_to_s[i][j][k]*vehicles[k][1] for k in range(len(vehicles))]) > 0:
                return False
    return True

if __name__ == "__main__":
    # 点的范围及个数
    x_min, y_min, x_max, y_max, nums_w, nums_d, nums_s = 0, 0, 100, 100, 5, 2, 5
    # 点集、运单生成
    points = points_generate(x_min, y_min, x_max, y_max, nums_w+nums_d+nums_s)
    ODs = OD_generate(1, 10000, nums_w, nums_d, nums_s)
    # points = [[0, 40, 39], [1, 99, 90], [2, 3, 85], [3, 91, 67], [4, 12, 36], [5, 31, 73], [6, 15, 22], [7, 64, 69]]
    # ODs = [[0, 0, 0, 0, 0, 2121, 1622, 379], [0, 0, 0, 0, 0, 7749, 877, 6880], [0, 0, 0, 0, 0, 1485, 1277, 700], [0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0]]

    print(points)
    print(ODs)
    OD_plot(points, ODs, x_min, y_min, x_max, y_max)

    # 约束
    constraint_w_to_d, constraint_d_to_s, capacity = constraint_generate(nums_w, nums_d, nums_s, nums_w*nums_s*10000/nums_s, nums_w*nums_s*10000*2/nums_s)
    # constraint_w_to_d = [[0, 2], [2, 0], [2, 0]]
    # constraint_d_to_s = [[0, 0, 1], [2, 1, 0]]
    # capacity = [80481, 85113]

    print(constraint_w_to_d)
    print(constraint_d_to_s)
    print(capacity)

    # 车辆（序号，载容，固定成本，变动成本）
    vehicles = [[0, 2000, 100, 4], [1, 1500, 60, 3]]

    # 生成路径
    R, ODs1 = R_generate_mix_cost_constraint(points, nums_w, nums_d, nums_s, ODs, vehicles, constraint_w_to_d, constraint_d_to_s, capacity, 100)
    OD_plot(points, ODs1, x_min, y_min, x_max, y_max)

    #3*2*2*3
    # R = [[[[0 for l in range(nums_s)] for k in range(nums_d)] for j in range(nums_d)] for i in range(nums_w)]
    # R[0][1][1][0] = 1
    # R[0][1][1][1] = 1
    # R[0][1][0][2] = 1
    # R[1][0][1][0] = 1
    # R[1][0][1][1] = 1
    # R[1][0][0][2] = 1
    # R[2][0][1][0] = 1
    # R[2][0][1][1] = 1
    # R[2][0][0][2] = 1
    #
    # bills_route_w_to_d = [[0,2121+1622+379],[7749+877+6880,0],[1485+1277+700,0]]
    # nums_vehicles_w_to_d = [[[0,0],[10,10]],[[10,10],[0,0]],[[10,10],[0,0]]]
    #
    # bills_route_d_to_d = [[6880+700, 7749+877+1485+1277], [379, 2121+1622]]
    # nums_vehicles_d_to_d = [[[0, 0], [10, 10]], [[10, 10], [0, 0]]]
    #
    # bills_route_d_to_s = [[0,0,379+6880+700], [2121+7749+1485,1622+877+1277,0]]
    # nums_vehicles_d_to_s = [[[0, 0], [0, 0], [10, 10]], [[10, 10], [10, 10], [0, 0]]]
    #
    # a = check_constraint(nums_w, nums_d, nums_s, bills_route_w_to_d, nums_vehicles_w_to_d, bills_route_d_to_d, nums_vehicles_d_to_d, bills_route_d_to_s, nums_vehicles_d_to_s,
    #                      constraint_w_to_d, constraint_d_to_s, capacity, R)
    #
    # print(a)