38、带约束的选址路径问题2

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)])
    LL_print("points", points)
    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+2)] for i in range(nums_w+nums_d+nums_s+2)]

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

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

def sub_LL_print(L, digit):
    s = "["
    for i in range(len(L)):
        s = s + str(L[i])
        temp = L[i]
        digit_temp = 1
        while temp//10>0:
            digit_temp += 1
            temp = temp//10
        for j in range(digit-digit_temp):
            s += " "
        if i < len(L)-1:
            s += ",  "
    s += "]"
    return s

def LL_print(name, LL):
    LL1 = [[int(LL[i][j]) for j in range(len(LL[i]))] for i in range(len(LL))]
    num_max = max([max(LL1[i]) for i in range(len(LL1))])
    digit = 1
    while num_max//10>0:
        digit += 1
        num_max = num_max // 10

    print("\n",name+":")
    print("["+sub_LL_print(LL[0], digit)+",")
    for i in range(1,len(LL)-1):
        print(" "+sub_LL_print(LL[i], digit)+",")
    print(" " + sub_LL_print(LL[-1], digit) + "]")

def constraint_w_to_d_generate(nums_w, nums_d):
    constraint_w_to_d = [[2 for j in range(nums_d)] for i in range(nums_w)]
    LL_print("constraint_w_to_d", constraint_w_to_d)
    return constraint_w_to_d

def constraint_d_to_s_generate(nums_d, nums_s):
    constraint_d_to_s = [[2 for j in range(nums_s)] for i in range(nums_d)]
    LL_print("constraint_d_to_s", constraint_d_to_s)
    return constraint_d_to_s

def capacity_generate(nums_d, capacity_low, capacity_up):
    capacity = [int((capacity_up - capacity_low) * rd.random() + capacity_up) for i in range(nums_d)]
    print("\ncapacity\n", capacity)
    return capacity

def in_in_generate(nums_d):
    in_in = [2 for i in range(nums_d)]
    print("\nin_in\n", in_in)
    return in_in

def in_out_generate(nums_d):
    in_out = [2 for i in range(nums_d)]
    print("\nin_out\n", in_out)
    return in_out

def out_in_generate(nums_d):
    out_in = [2 for i in range(nums_d)]
    print("\nout_in\n", out_in)
    return out_in

def out_out_generate(nums_d):
    out_out = [2 for i in range(nums_d)]
    print("\nout_out\n", out_out)
    return out_out

def rent_per_generate(nums_d):
    rent_per = [0.68 for i in range(nums_d)]
    print("\nrent_per\n", rent_per)
    return rent_per

def nums_limit_generate(nums_d):
    nums_limit = [nums_d, int(nums_d/2), nums_d, int(nums_d/2)]
    print("\nnums_limit\n", nums_limit)
    return nums_limit

def vehicles_generate(name):
    vehicles = [[0, 2000, 100, 4], [1, 1500, 60, 3]]
    LL_print(name, vehicles)
    return vehicles

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_w_to_d, vehicles_d_to_d, vehicles_d_to_s, constraint_w_to_d, constraint_d_to_s,
                                   capacity, in_in, in_out, out_in, out_out,
                                   rent_per, nums_limit, time_limit, OD_sum):
    # 欧式距离和球面距离
    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)]
    LL_print("Dis", Dis)

    model = Model("Level_mix_cost_constraint_2")

    M = sum([sum(ODs[i]) for i in range(len(ODs))])*10
    print("\nM:\n", M)

    R = [[[[model.addVar(vtype="B", name="R[%s,%s,%s,%s]" % (i, j, k, l)) for l in range(nums_s+2)] for k in range(nums_d)] for j in range(nums_d)] for i in range(nums_w+nums_s+2)]

    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+nums_s+2)]
    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+2)] 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_w_to_d))] for j in range(nums_d)] for i in range(nums_w+nums_s+2)]
    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_d_to_d))] 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_d_to_s))] for j in range(nums_s+2)] for i in range(nums_d)]

    bills_d = [model.addVar(vtype="I", name="bills_d[%s]" % i) for i in range(nums_d)]
    flag_d = [[model.addVar(vtype="B", name="flag_d[%s,%s]" % (i, j)) for j in range(4)] for i in range(nums_d)]

    # 以总成本最小为目标
    model.setObjective(quicksum(bills_d[i]*rent_per[i] for i in range(nums_d)) + quicksum((Dis[i][j+nums_w] * vehicles_w_to_d[k][3] + vehicles_w_to_d[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_w_to_d))) + quicksum((Dis[i+nums_w][j+nums_w] * vehicles_d_to_d[k][3] + vehicles_d_to_d[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_d_to_d)) if i != j) + quicksum((Dis[i+nums_w][j+nums_w+nums_d] * vehicles_d_to_s[k][3] + vehicles_d_to_s[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_d_to_s))), "minimize")

    # 变量约束
    for i in range(nums_w+nums_s+2):
        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_w_to_d)):
                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_d_to_d)):
                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+2):
            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_d_to_s)):
                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_d):
        model.addCons(bills_d[i] >= 0)
        model.addCons(bills_d[i] <= M)

    for i in range(nums_d):
        for j in range(4):
            model.addCons(flag_d[i][j] >= 0)
            model.addCons(flag_d[i][j] <= 1)

    #单条路线约束
    for i in range(nums_w+nums_s+2):
        for l in range(nums_s+2):
            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 i in range(nums_s):
        flag = 0
        idex_delv = 0
        for j in range(nums_d):
            if constraint_d_to_s[j][i] == 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+nums_w][j] == 0)
        else:
            for j in range(nums_d):
                if constraint_d_to_s[j][i] == 0:
                    model.addCons(bills_route_w_to_d[i+nums_w][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) if i != j) + quicksum(bills_route_d_to_s[i][j] for j in range(nums_s)) <= capacity[i])

    #运输单量平衡
    for i in range(nums_w+nums_s+2):
        for j in range(nums_d):
            if i < nums_w:
                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+2)) == 0)
            else:
                model.addCons(bills_route_w_to_d[i][j] - quicksum(ODs[i+nums_d][l + nums_w + nums_d] * R[i][j][k][l] for k in range(nums_d) for l in range(nums_s + 2)) == 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+2))-quicksum(ODs[i+nums_d][l+nums_w+nums_d]*R[i][j][k][l] for i in range(nums_w, nums_w+nums_s+2) for l in range(nums_s+2)) == 0)

    for k in range(nums_d):
        for l in range(nums_s+2):
            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))-quicksum(ODs[i+nums_d][l+nums_w+nums_d]*R[i][j][k][l] for i in range(nums_w, nums_w+nums_s+2) 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_w_to_d[k][1] for k in range(len(vehicles_w_to_d))) <= 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_d_to_d[k][1] for k in range(len(vehicles_d_to_d))) <= 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_d_to_s[k][1] for k in range(len(vehicles_d_to_s))) <= 0)


    #分拣能力约束
    for j in range(nums_d):
        if in_in[j] == 0:
            model.addCons(quicksum(R[nums_w+nums_s][j][k][l] for k in range(nums_d) for l in range(nums_s+2)) == 0)
        if in_out[j] == 0:
            model.addCons(quicksum(R[nums_w+nums_s+1][j][k][l] for k in range(nums_d) for l in range(nums_s+2)) == 0)

    for k in range(nums_d):
        if out_in[k] == 0:
            model.addCons(quicksum(R[i][j][k][nums_w + nums_s] for i in range(nums_w + nums_s + 2) for l in range(nums_d)) == 0)
        if out_out[k] == 0:
            model.addCons(quicksum(R[i][j][k][nums_w + nums_s + 1] for i in range(nums_w + nums_s + 2) for l in range(nums_d)) == 0)


    #运单与示性函数约束
    for i in range(nums_d):
        model.addCons(bills_route_w_to_d[nums_w + nums_s][i] - M * flag_d[i][0] <= 0)
        model.addCons(bills_route_w_to_d[nums_w + nums_s + 1][i] - M * flag_d[i][1] <= 0)
        model.addCons(bills_route_d_to_s[i][nums_s] - M * flag_d[i][2] <= 0)
        model.addCons(bills_route_d_to_s[i][nums_s + 1] - M * flag_d[i][3] <= 0)

    #分拣功能个数约束
    model.addCons(quicksum(flag_d[i][0] for i in range(nums_d)) - nums_limit[0] <= 0)
    model.addCons(quicksum(flag_d[i][1] for i in range(nums_d)) - nums_limit[1] <= 0)
    model.addCons(quicksum(flag_d[i][2] for i in range(nums_d)) - nums_limit[2] <= 0)
    model.addCons(quicksum(flag_d[i][3] for i in range(nums_d)) - nums_limit[3] <= 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+2)] for k in range(nums_d)] for j in range(nums_d)] for i in range(nums_w+nums_s+2)]
    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+nums_s+2)]
    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+2)] for i in range(nums_d)]

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

    #车辆运载约束
    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+2):
            ODs1[i+nums_w][j+nums_w+nums_d] = bills_route_d_to_s1[i][j]

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

    #打印结果路线
    print("\n路线")
    for i in range(nums_w+nums_s+2):
        for j in range(nums_d):
            for k in range(nums_d):
                for l in range(nums_s+2):
                    if R1[i][j][k][l] == 1:
                        print(i,"-",j,"-",k,"-",l)

    LL_print("bills_route_w_to_d1", bills_route_w_to_d1)
    LL_print("bills_route_d_to_d1", bills_route_d_to_d1)
    LL_print("bills_route_d_to_s1", bills_route_d_to_s1)

    return R1, ODs1, bills_route_d_to_d1

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()

if __name__ == "__main__":
    # 点的范围及个数
    x_min, y_min, x_max, y_max, nums_w, nums_d, nums_s = 0, 0, 100, 100, 1, 2, 1
    # 点集、运单生成
    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)
    OD_plot(points, ODs, x_min, y_min, x_max, y_max)
    # 总单量
    OD_sum = sum([sum(ODs[i]) for i in range(nums_w)])
    print("\n总单量:\n" , OD_sum)
    # 添加约束（仓拣、拣站、产能、支线进港、干线进港、支线出港、干线出港、站点揽收、站点派送、操作成本、个数上限）
    constraint_w_to_d = constraint_w_to_d_generate(nums_w, nums_d)
    constraint_d_to_s = constraint_d_to_s_generate(nums_d, nums_s)
    capacity = capacity_generate(nums_d, OD_sum, OD_sum * 2)
    in_in = in_in_generate(nums_d)
    in_out = in_out_generate(nums_d)
    out_in = out_in_generate(nums_d)
    out_out = out_out_generate(nums_d)
    rent_per = rent_per_generate(nums_d)
    nums_limit = nums_limit_generate(nums_d)
    # 车辆
    vehicles_w_to_d = vehicles_generate("vehicles_w_to_d")
    vehicles_d_to_d = vehicles_generate("vehicles_d_to_d")
    vehicles_d_to_s = vehicles_generate("vehicles_d_to_s")

    # 数据准备好，开始计算
    R, ODs1, bills_route_d_to_d1 = R_generate_mix_cost_constraint(points, nums_w, nums_d, nums_s, ODs, vehicles_w_to_d, vehicles_d_to_d, vehicles_d_to_s,
                                       constraint_w_to_d, constraint_d_to_s, capacity, in_in, in_out, out_in, out_out,
                                       rent_per, nums_limit, 1000, OD_sum)
    OD_plot(points, ODs1, x_min, y_min, x_max, y_max)