30、二维装箱（单品）

from pyscipopt import Model, quicksum
from vtk import *
import vtk
import random as rd
import time
import numpy as np
import functools
import copy

#数据生成，输入为箱子种类数，箱子最大长、最小长、最大宽、最小宽、商品个数、商品最大长、最小长、最大宽、最小宽
#输出为生成的箱子列表和商品列表
def data_generate_2d_simple(nums_box, max_edge_box, min_edge_box, nums_good, max_edge_good, min_edge_good):

    #随机生成箱子
    boxs = []
    for i in range(nums_box):
        boxs.append([int((max_edge_box-min_edge_box) * rd.random() + min_edge_box), int((max_edge_box-min_edge_box) * rd.random() + min_edge_box)])

    #随机生成商品
    goods = []
    good = [int((max_edge_good-min_edge_good) * rd.random() + min_edge_good), int((max_edge_good-min_edge_good) * rd.random() + min_edge_good)]
    for i in range(nums_good):
        goods.append(copy.deepcopy(good))

    return boxs,goods

#二维商品排序用
def cmp_2d(x,y):
    if x[0] < y[0]:
        return -1
    elif x[0] > y[0]:
        return 1
    elif x[1] < y[1]:
        return -1
    elif x[1] > y[1]:
        return 1
    else:
        return 0

def data_pre(boxs, goods):

    #箱子长边作长，放前面
    for i in range(len(boxs)):
        boxs[i] = sorted(boxs[i], reverse=True)

    #箱子排序，大的放前面
    boxs = sorted(boxs, key=functools.cmp_to_key(cmp_2d), reverse=True)

    #箱子种类去重
    i = 0
    while i < len(boxs) - 1:
        j = i + 1
        while j < len(boxs):
            if boxs[i][0] == boxs[j][0] and boxs[i][1] == boxs[j][1]:
                del boxs[j]
            else:
                j += 1
        i += 1

    #商品长边作长，放前面
    for i in range(len(goods)):
        goods[i] = sorted(goods[i], reverse=True)

    #商品排序，大的放前面
    goods = sorted(goods, key=functools.cmp_to_key(cmp_2d), reverse=True)

    return boxs,goods

#检验是否每个商品都能有一个箱子放下它
def check_2d(boxs, goods):
    for good in goods:
        can_put = False
        for box in boxs:
            if good[0] <= box[0] and good[1] <= box[1]:
                can_put = True
                break
        if not can_put:
            print(good,"太大，无合适箱子")
            return False
    return True

#单品类五块方案
def five_block_simple_2d(CL, CW, l, w):

    maxnum = int(CL/l) * int(CW/w) + int((CL - int(CL/l)*l)/w) * int(CW/l)
    sub_L_coordinate = [[l,w,l*i,w*j] for i in range(int(CL/l)) for j in range(int(CW/w))] + \
                        [[w, l, int(CL/l)*l + w * i, l * j] for i in range(int((CL - int(CL/l)*l)/w)) for j in range(int(CW/l))]
    maxnum2 = int(CW/l) * int(CL/w) + int((CW - int(CW/l)*l)/w) * int(CL/l)
    sub_L_coordinate2 = [[w, l, w * i, l * j] for i in range(int(CL / w)) for j in range(int(CW / l))] + \
                         [[l, w, l * i, int(CW/l)*l + w * j] for i in range(int(CL/l)) for j in range(int((CW - int(CW/l)*l)/w))]

    if maxnum < maxnum2:
        maxnum, sub_L_coordinate = maxnum2, sub_L_coordinate2

    #五块调整, 当剩余面积大于一块底面积时触发
    if int(CL * CW / l / w) > maxnum and CW >= l + w:
        up1, up2, up3, up4 = int(CL / l), int(CL / w), int(CW / l), int(CW / w)
        for x1 in range(1, up1 + 1):
            x2 = int((CL - x1 * l) / w)
            for y1 in range(1, up4 + 1):
                y4 = int((CW-y1 * w) / l)
                for y2 in range(1, up3 + 1):
                    y3 = int((CW-y2 * l) / w)
                    for x3 in range(1, up1 + 1):
                        x4 = int((CL-x3 * l) / w)
                        if y1 * w <= y2 * l and x2 * w <= x3 * l and y3 * w <= y4 * l and x4 * w <= x1 * l:
                            if ((y1 + y3) * w > CW or (x2 + x4) * w > CL):
                                continue
                            n = x1 * y1 + x2 * y2 + x3 * y3 + x4 * y4
                            temp_sub_L_coordinate = [[l, w, l * i, w * j] for i in range(x1) for j in range(y1)] + \
                                                     [[w, l, CL-w*(i+1), l*j] for i in range(x2) for j in range(y2)] + \
                                                     [[l, w, CL - l*(i+1), CW - w*(j+1)] for i in range(x3) for j in range(y3)] + \
                                                     [[w, l, w*i, CW - l*(j+1)] for i in range(x4) for j in range(y4)]
                            if maxnum < n:
                                maxnum, sub_L_coordinate = n, temp_sub_L_coordinate

                        if y1 * w >= y2 * l and x2 * w >= x3 * l and y3 * w >= y4 * l and x4 * w >= x1 * l:
                            if (x1+x3) * l > CL or (y2+y4) * l > CW:
                                continue
                            n = x1 * y1 + x2 * y2 + x3 * y3 + x4 * y4
                            temp_sub_L_coordinate = [[l, w, l * i, w * j] for i in range(x1) for j in range(y1)] + \
                                                     [[w, l, CL - w * (i + 1), l * j] for i in range(x2) for j in range(y2)] + \
                                                     [[l, w, CL - l * (i + 1), CW - w * (j + 1)] for i in range(x3) for j in range(y3)] + \
                                                     [[w, l, w * i, CW - l * (j + 1)] for i in range(x4) for j in range(y4)]
                            if maxnum < n:
                                maxnum, sub_L_coordinate = n, temp_sub_L_coordinate

    return [maxnum, sub_L_coordinate]

#箱子选择，V、W：每个箱子的容纳个数、体积，b：总个数
def bin_choose(V, W, b, time_limit):

    n = len(V)
    model = Model("bin_choose")
    #每个箱子选几个
    X = [model.addVar(vtype="I", name="x[%s]" % i) for i in range(n)]

    #以总包裹体积最小（即填充率最大）为目标
    model.setObjective(quicksum(W[i]*X[i] for i in range(n)), "minimize")

    #装载个数大于b
    model.addCons(quicksum(V[i]*X[i] for i in range(n)) - b >= 0)

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

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

    X1 = [round(model.getVal(X[i])) for i in range(n)]

    return X1

#多品类五块方案
def five_block_multiple_2d(boxs, goods, time_limit):

    #分别计算每个箱子能装几个商品
    nums, sub_L_coordinates = [], []
    for i in range(len(boxs)):
        r = five_block_simple_2d(boxs[i][0], boxs[i][1], goods[0][0], goods[0][1])
        nums.append(r[0])
        sub_L_coordinates.append(r[1])

    X = bin_choose(nums, [boxs[i][0]*boxs[i][1] for i in range(len(boxs))], len(goods), time_limit)

    L_box = [boxs[i] for i in range(len(boxs)) for j in range(X[i])]
    L_goods = [[goods[0] for k in range(nums[i])] for i in range(len(boxs)) for j in range(X[i])]
    L_coordinates = [sub_L_coordinates[i] for i in range(len(boxs)) for j in range(X[i])]
    over = sum([X[i]*nums[i] for i in range(len(X))]) - len(goods)
    for i in range(len(L_goods)):
        if len(L_goods[i]) >= over:
            L_goods[i] = L_goods[i][:len(L_goods[i])-over]
            L_coordinates[i] = copy.deepcopy(L_coordinates[i][:len(L_coordinates[i]) - over])
            break

    return L_box, L_goods, L_coordinates


#检验结果中的商品集是否和原始的商品集一致
def goods_check(goods, L_goods):

    nums = 0
    for gs in L_goods:
        nums += len(gs)

    if len(goods) == nums:
        return True

    return False

#二维装箱整数规划
def IP_2d(boxs, goods, time_limit):

    m = len(boxs)
    n = len(goods)

    M = max([max(boxs[i]) for i in range(m)])*10

    model = Model("IP_2d")

    #对X、Y，i为包裹下标，j为商品下标，x[i][j]=1代表第j件商品放进第i个包裹里，y[i]=1表示开启第i个包裹
    X = [[model.addVar(vtype="B", name="X[%s,%s]" % (i, j)) for j in range(n)] for i in range(m*n)]
    Y = [model.addVar(vtype="B", name="Y[%s]" % i) for i in range(m*n)]
    #包裹面积向量，[n个boxs[0],n个boxs[1],...,n个boxs[m-1]]
    c = [boxs[i][0]*boxs[i][1] for i in range(len(boxs)) for j in range(len(goods))]
    cb = [boxs[i] for i in range(len(boxs)) for j in range(len(goods))]

    #码放形式，1代表(l,w), 0代表(w,l)
    type_place = [model.addVar(vtype="B", name="tp[%s]" % i) for i in range(n)]

    #码放长、宽、坐标x、y、z
    L = [model.addVar(vtype="I", name="L[%s]" % i) for i in range(n)]
    W = [model.addVar(vtype="I", name="W[%s]" % i) for i in range(n)]
    x = [model.addVar(vtype="I", name="x[%s]" % i) for i in range(n)]
    y = [model.addVar(vtype="I", name="y[%s]" % i) for i in range(n)]

    #选择变量，fx[i][j]=1代表i在j的左边（小），fy[i][j]=1同理，f[i][j][k]=1代表j、k同时在第i个箱子中
    fx = [[model.addVar(vtype="B", name="fx[%s,%s]" % (i, j)) for j in range(n)] for i in range(n)]
    fy = [[model.addVar(vtype="B", name="fy[%s,%s]" % (i, j)) for j in range(n)] for i in range(n)]
    f = [[[model.addVar(vtype="B", name="f[%s,%s,%s]" % (i, j, k)) for k in range(n)] for j in range(n)] for i in range(m*n)]

    #以总包裹体积最小（即填充率最大）为目标
    model.setObjective(quicksum(Y[i]*c[i] for i in range(m*n)), "minimize")

    # 使用前必须开启
    for i in range(m * n):
        for j in range(n):
            model.addCons(X[i][j] - Y[i] <= 0)

    #每个商品都能被装下
    for j in range(n):
        model.addCons(quicksum(X[i][j] for i in range(m * n)) == 1)

    #码放长宽约束
    for j in range(n):
        model.addCons(L[j] - goods[j][0] * type_place[j] - goods[j][1] * (1 - type_place[j]) == 0)
        model.addCons(W[j] - goods[j][1] * type_place[j] - goods[j][0] * (1 - type_place[j]) == 0)

    #位置约束，大于0，不超过边界，不相交
    for j in range(n):
        model.addCons(x[j] >= 0)
        model.addCons(y[j] >= 0)

    for i in range(m*n):
        for j in range(n):
            model.addCons(x[j] + L[j] - cb[i][0] - M*(1 - X[i][j]) <= 0)
            model.addCons(y[j] + W[j] - cb[i][1] - M*(1 - X[i][j]) <= 0)

    for j in range(n):
        for k in range(j+1,n):
            model.addCons(x[j] + L[j] - x[k] - M*(1 - fx[j][k]) <= 0)
            model.addCons(x[k] + L[k] - x[j] - M*(1 - fx[k][j]) <= 0)
            model.addCons(y[j] + W[j] - y[k] - M*(1 - fy[j][k]) <= 0)
            model.addCons(y[k] + W[k] - y[j] - M*(1 - fy[k][j]) <= 0)

    for i in range(m*n):
        for j in range(n):
            for k in range(j+1,n):
                model.addCons(X[i][j] + X[i][k] - 1 - M*(1 - f[i][j][k]) <= 0)

    for j in range(n):
        for k in range(j+1,n):
            model.addCons(fx[j][k] + fx[k][j] + fy[j][k] + fy[k][j] + quicksum(f[i][j][k] for i in range(m*n)) >= m*n)

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

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

    #拿结果
    X1 = [[round(model.getVal(X[i][j])) for j in range(n)] for i in range(m * n)]
    L1 = [round(model.getVal(L[i])) for i in range(n)]
    W1 = [round(model.getVal(W[i])) for i in range(n)]
    x1 = [round(model.getVal(x[i])) for i in range(n)]
    y1 = [round(model.getVal(y[i])) for i in range(n)]

    L_box = []
    L_goods = []
    L_coordinates = []

    for i in range(m*n):
        goods_i = []
        coordinates_i = []
        for j in range(n):
            if X1[i][j] == 1:
                goods_i.append(goods[j])
                coordinates_i.append([L1[j],W1[j],x1[j],y1[j]])
        if len(goods_i) > 0:
            L_box.append(cb[i])
            L_goods.append(goods_i)
            L_coordinates.append(coordinates_i)

    return L_box, L_goods, L_coordinates, model.getGap()

#任务分流汇总
#给一系列箱子和商品（箱子可用个数不限），推荐结果
def stacking_2d_simple(boxs, goods, time_limit, nums_limit):

    #长宽预处理，降序排序，箱子去重
    boxs, goods = data_pre(boxs, goods)

    #这里是否所有的商品均至少有一个箱子可以装下，若有商品超出规格则直接返回
    if not check_2d(boxs, goods):
        return [[], []]

    #当商品数超过一定数量时，直接采用启发式算法
    if len(goods) > nums_limit:
        return five_block_multiple_2d(boxs, goods, time_limit)

    # 多类箱子，应用整数规划求解
    L_box, L_goods, L_coordinates, gap = IP_2d(boxs, goods, time_limit)

    # 结果检验，当求解器的结果有问题（商品数不符时）采用混合OBT求解方案
    if not goods_check(goods, L_goods):
        return five_block_multiple_2d(boxs, goods, time_limit)

    # gap较大时，用启发式方法比较，取优，相同时采用启发式方法（较为规整）
    if gap >= 0.01:
        print("尝试采用启发式方法")
        L_box1, L_goods1, L_coordinates1 = five_block_multiple_2d(boxs, goods, time_limit)
        if sum([L_box[i][0]*L_box[i][1] for i in range(len(L_box))]) >= sum([L_box1[i][0]*L_box1[i][1] for i in range(len(L_box1))]) and goods_check(goods, L_goods1):
            print("采用启发式方法")
            L_box, L_goods, L_coordinates = L_box1, L_goods1, L_coordinates1

    return L_box, L_goods, L_coordinates

#添加商品图形
def Addcube_2d_simple(ren, coordinate, edge_max, h, x_re, y_re, z_re):
    cube = vtk.vtkCubeSource()
    cube.SetXLength(coordinate[0]/edge_max)
    cube.SetYLength(coordinate[1]/edge_max)
    cube.SetZLength(h)
    cube.Update()

    translation = vtkTransform()
    translation.Translate((coordinate[2] + coordinate[0]/2.0)/edge_max + x_re, (coordinate[3] + coordinate[1]/2.0)/edge_max + y_re, h/2.0 + z_re)
    transformFilter = vtkTransformPolyDataFilter()
    transformFilter.SetInputConnection(cube.GetOutputPort())
    transformFilter.SetTransform(translation)
    transformFilter.Update()

    transformedMapper = vtkPolyDataMapper()
    transformedMapper.SetInputConnection(transformFilter.GetOutputPort())
    transformedActor = vtkActor()
    transformedActor.SetMapper(transformedMapper)
    transformedActor.GetProperty().SetColor((rd.uniform(0, 1), rd.uniform(0, 1), rd.uniform(0, 1)))

    ren.AddActor(transformedActor)

def png_save(renWin, name):
    windowToImageFilter = vtkWindowToImageFilter()
    windowToImageFilter.SetInput(renWin)
    windowToImageFilter.Update()
    writer = vtkPNGWriter()
    writer.SetFileName(name)
    writer.SetInputConnection(windowToImageFilter.GetOutputPort())
    writer.Write()

#二维展示，输入为箱子集和商品集，包裹的箱子和商品集一一对应
def show_2d_simple(L_box, L_coordinates):

    nums = len(L_box)
    edge_max = max([max(L_box[i]) for i in range(len(L_box))]) if max([max(L_box[i]) for i in range(len(L_box))]) > 0 else 1

    #预设参数
    CH = 0.5
    gap = 0.25
    CL_p = 1.1
    CW_p = nums + gap * (nums - 1)
    CH_p = 0.01
    gap = 0.25

    x_re = -0.5
    y_re = -0.5
    z_re = -0.5

    #渲染及渲染窗口，并根据捕捉的鼠标事件执行相应的操作
    ren = vtk.vtkRenderer()
    renWin = vtk.vtkRenderWindow()
    renWin.AddRenderer(ren)
    renWin.SetSize(1200, 600)
    iren = vtk.vtkRenderWindowInteractor()
    iren.SetRenderWindow(renWin)


    """画容器"""
    for i in range(nums):

        cube = vtk.vtkCubeSource()
        cube.SetXLength(L_box[i][0]/edge_max)
        cube.SetYLength(L_box[i][1]/edge_max)
        cube.SetZLength(CH)
        cube.Update()

        translation = vtkTransform()
        translation.Translate(L_box[i][0]/edge_max/2.0 + x_re, L_box[i][1]/edge_max/2.0 + i + gap*i + y_re, CH/2.0 + z_re)
        transformFilter = vtkTransformPolyDataFilter()
        transformFilter.SetInputConnection(cube.GetOutputPort())
        transformFilter.SetTransform(translation)
        transformFilter.Update()

        transformedMapper = vtkPolyDataMapper()
        transformedMapper.SetInputConnection(transformFilter.GetOutputPort())
        transformedActor = vtkActor()
        transformedActor.SetMapper(transformedMapper)
        transformedActor.GetProperty().SetColor((1, 1, 1))
        transformedActor.GetProperty().SetRepresentationToWireframe()

        ren.AddActor(transformedActor)

    """画托盘"""
    cube = vtk.vtkCubeSource()
    cube.SetXLength(CL_p)
    cube.SetYLength(CW_p)
    cube.SetZLength(CH_p)
    cube.Update()

    translation = vtkTransform()
    translation.Translate(CL_p/2.0 + x_re, CW_p/2.0 + y_re, -CH_p/2.0 + z_re)
    transformFilter = vtkTransformPolyDataFilter()
    transformFilter.SetInputConnection(cube.GetOutputPort())
    transformFilter.SetTransform(translation)
    transformFilter.Update()

    transformedMapper = vtkPolyDataMapper()
    transformedMapper.SetInputConnection(transformFilter.GetOutputPort())
    transformedActor = vtkActor()
    transformedActor.SetMapper(transformedMapper)
    transformedActor.GetProperty().SetColor((0.2, 0.4, 0.8))

    ren.AddActor(transformedActor)

    for i in range(len(L_coordinates)):
        for j in range(len(L_coordinates[i])):
            Addcube_2d_simple(ren, L_coordinates[i][j], edge_max, CH, x_re, i + gap*i + y_re, z_re)

    camera = vtk.vtkCamera()
    camera.SetPosition(5, -0.5, 2)
    camera.SetViewUp(0, 0, 1)
    ren.SetActiveCamera(camera)

    iren.Initialize()
    renWin.Render()
    # 保存过程
    png_save(renWin, "result_D2_simple.png")
    # 展示
    iren.Start()

# 数值实验， 输入为箱子种类数，箱子最大容积、最小容积、商品个数、商品最大体积、商品最小体积、时间限制、实验次数，
# 输出为five_block平均装载率、OBT平均时间、整数规划平均装载率、整数规划平均时间
def experiment_2d_simple(nums_box, max_edge_box, min_edge_box, nums_good, max_edge_good, min_edge_good, time_limit, times_experiment):

    rate_five_block = []
    time_five_block = []
    rate_IP = []
    time_IP = []

    for i in range(times_experiment):
        for j in range(20):
            print("--------------------------", nums_good, "----------------------", i)
        boxs, goods = data_generate_2d_simple(nums_box, max_edge_box, min_edge_box, nums_good, max_edge_good, min_edge_good)
        # 长宽预处理，降序排序，箱子去重
        boxs, goods = data_pre(boxs, goods)
        t1 = time.clock()
        L_box1, L_goods1, L_coordinates1 = five_block_multiple_2d(boxs, goods, time_limit)
        t2 = time.clock()
        L_box2, L_goods2, L_coordinates2, gap = IP_2d(boxs, goods, time_limit)
        t3 = time.clock()
        if sum([L_box1[i][0]*L_box1[i][1] for i in range(len(L_box1))]) > 0 and goods_check(goods, L_goods1):
            rate_five_block.append(sum([goods[i][0]*goods[i][1] for i in range(len(goods))])/sum([L_box1[i][0]*L_box1[i][1] for i in range(len(L_box1))]))
        else:
            rate_five_block.append(0)

        time_five_block.append(t2-t1)
        if sum([L_box2[i][0]*L_box2[i][1] for i in range(len(L_box2))]) > 0 and goods_check(goods, L_goods2):
            rate_IP.append(sum([goods[i][0]*goods[i][1] for i in range(len(goods))])/sum([L_box2[i][0]*L_box2[i][1] for i in range(len(L_box2))]))
        else:
            rate_IP.append(0)
        time_IP.append(t3-t2)

    print("rate_five_block:", np.mean(rate_five_block))
    print("rate_IP:", np.mean(rate_IP))
    print("time_five_block:", np.mean(time_five_block))
    print("time_IP:", np.mean(time_IP))

    return np.mean(rate_five_block), np.mean(rate_IP), np.mean(time_five_block), np.mean(time_IP)

if __name__ == "__main__":

    # # 生成箱子集和商品集，计算并展示
    # boxs,goods = data_generate_2d_simple(nums_box = 2, max_edge_box = 20, min_edge_box = 10, nums_good = 10, max_edge_good = 10, min_edge_good = 1)
    # # boxs, goods = [[19, 18], [15, 12]], [[8, 6], [8, 6], [8, 6], [8, 6], [8, 6], [8, 6], [8, 6], [8, 6], [8, 6], [8, 6]]
    # print(boxs)
    # print(goods)
    # L_box, L_goods, L_coordinates = stacking_2d_simple(boxs, goods, time_limit=100, nums_limit=50)
    # print(L_box)
    # print(L_coordinates)
    # show_2d_simple(L_box, L_coordinates)

    # 数值实验
    rate_five_block = []
    rate_IP = []
    time_five_block = []
    time_IP = []
    for i in range(30):
        r1, r2, r3, r4 = experiment_2d_simple(nums_box=2, max_edge_box=20, min_edge_box=10, nums_good=i + 1, max_edge_good=10,
                                       min_edge_good=1, time_limit=100, times_experiment=50)
        rate_five_block.append(r1)
        rate_IP.append(r2)
        time_five_block.append(r3)
        time_IP.append(r4)
        print("-------------------------", i + 1, "个商品测试完成")
        print(rate_five_block)
        print(rate_IP)
        print(time_five_block)
        print(time_IP)