Advanced constraint programming modeling

advancedModelingWeek4.py

from ortools.sat.python import cp_model

''' course: https://www.coursera.org/learn/advanced-modeling#about '''



def squarePackingProblem():
    model = cp_model.CpModel()

    copyList = [4, 3, 0, 5, 4, 3, 4]
    numOfSquareSizes = len(copyList)
    maxLength = sum([(i+1) * copyList[i] for i in range(numOfSquareSizes)])
    minArea = sum([(i+1) * (i+1) * copyList[i] for i in range(numOfSquareSizes)])

    sizeList, firstIndexList, numOfSquares = [], [], 0
    for square in range(numOfSquareSizes):
        for _ in range(copyList[square]):
            sizeList.append(square+1)
            firstIndexList.append(numOfSquares)
        numOfSquares += copyList[square]


    heightVar = model.NewIntVar(numOfSquareSizes, maxLength, "height")
    widthVar = model.NewIntVar(numOfSquareSizes, maxLength, "width")
    areaVar = model.NewIntVar(minArea, numOfSquareSizes*maxLength, "area")

    model.AddMultiplicationEquality(areaVar, [widthVar, heightVar])

    xVarList = []
    xEndVarList = []
    yVarList = []
    yEndVarList = []
    for square in range(numOfSquares):
        xVarList.append(model.NewIntVar(0, maxLength, "xPositionForSquare"+str(square)))
        xEndVarList.append(model.NewIntVar(0, maxLength, "xEndForSquare" + str(square)))
        yVarList.append(model.NewIntVar(0, maxLength, "yPositionForSquare" + str(square)))
        yEndVarList.append(model.NewIntVar(0, maxLength, "yEndForSquare" + str(square)))

    for square in range(numOfSquares):
        model.Add(xVarList[square] + sizeList[square] <= widthVar)
        model.Add(yVarList[square] + sizeList[square] <= heightVar)

    xIntervalVarList = []
    yIntervalVarList = []
    for square in range(numOfSquares):
        xIntervalVarList.append(model.NewIntervalVar(xVarList[square], sizeList[square], xEndVarList[square], "xIntervalVarFor"+str(square)))
        yIntervalVarList.append(model.NewIntervalVar(yVarList[square], sizeList[square], yEndVarList[square], "yIntervalVarFor"+str(square)))


    model.AddCumulative(xIntervalVarList, sizeList, heightVar) # redundant but improves propagation
    model.AddCumulative(yIntervalVarList, sizeList, widthVar) # redundant but improves propagation
    model.AddNoOverlap2D(xIntervalVarList, yIntervalVarList)


    square = 0
    for _ in range(numOfSquareSizes):  # symmetry breaking using lexicographic ordering
        if copyList[square] <= 0: continue
        for copy in range(copyList[square] - 1):
            boolVar = model.NewBoolVar("boolVar"+str(square)+str(copy))
            model.Add(xVarList[firstIndexList[square] + copy] > xVarList[firstIndexList[square] + copy + 1]).OnlyEnforceIf(boolVar)
            model.Add(xVarList[firstIndexList[square] + copy] == xVarList[firstIndexList[square] + copy + 1]).OnlyEnforceIf(boolVar.Not())
            model.Add(yVarList[firstIndexList[square] + copy] > yVarList[firstIndexList[square] + copy + 1]).OnlyEnforceIf(boolVar.Not())
        square += 1


    xSol = [20, 20, 5, 5, 9, 4, 4, 21, 21, 5, 5, 0, 20, 6, 0, 0, 15, 9, 0, 18, 13, 11, 6]
    ySol = [6, 5, 15, 6, 13, 18, 16, 9, 5, 11, 7, 16, 0, 15, 11, 6, 7, 7, 0, 13, 0, 13, 0]
    for i,var1 in enumerate(xVarList):
        model.AddHint(var1, xSol[i])
    for i,var2 in enumerate(yVarList):
        model.AddHint(var2, ySol[i])

    model.Minimize(areaVar)

    solver = cp_model.CpSolver()
    status = solver.Solve(model)

    if status == cp_model.INFEASIBLE:
        print("Infeasible model")

    if status == cp_model.OPTIMAL:
        print('Minimum area: %i' % solver.ObjectiveValue())
        print("height:", solver.Value(heightVar))
        print("width:", solver.Value(widthVar))
        print("x =", [solver.Value(var) for var in xVarList])
        print("y =", [solver.Value(var) for var in yVarList])
    print("----------------------------")




#----------------------------------------------------------------------------------------------------------------------#




def rectilinearPackingProblem():
    model = cp_model.CpModel()

    offsetList = [(1,0,2,5), (3,4,1,1), (0,3,1,1), (1,4,2,2), (0,1,1,3), (1,0,4,4), # xoffset,yoffset,xsize,ysize
                  (1,5,1,2), (1,0,3,5), (4,1,1,4), (0,0,1,3), (1,0,3,4), (0,0,4,5)]
    shapeList = [(0,1,2), (3,4,5), (4,6,7,8), (11,), (9,10)]
    numOfBlocks = len(shapeList)
    widthOfRiver = 9
    maxLength = 16


    xVarList = []
    yVarList = []
    for block in range(numOfBlocks):
        xVarList.append(model.NewIntVar(0, maxLength, "xPositionForBlock"+str(block)))
        yVarList.append(model.NewIntVar(0, widthOfRiver, "yPositionForBlock" + str(block)))

    lengthUsedVar = model.NewIntVar(0, maxLength, "lengthUsed")

    for block in range(numOfBlocks):
        for shape in shapeList[block]:
            offset = offsetList[shape]
            model.Add(xVarList[block] + offset[0] + offset[2] <= lengthUsedVar)
            model.Add(yVarList[block] + offset[1] + offset[3] <= widthOfRiver)


    for block1 in range(numOfBlocks-1):
        for shape1 in shapeList[block1]:
            offset1 = offsetList[shape1]

            for block2 in range(block1+1, numOfBlocks):
                for shape2 in shapeList[block2]:
                    offset2 = offsetList[shape2]

                    block1IsLeftOfBlock2 = model.NewBoolVar("blockShape"+str(block1)+str(shape1)+"isLeftTo"+str(block2)+str(shape2))
                    block2IsLeftOfBlock1 = model.NewBoolVar("blockShape"+str(block2)+str(shape2)+"isLeftTo"+str(block1)+str(shape1))
                    block1IsBelowBlock2 = model.NewBoolVar("blockShape"+str(block1)+str(shape1)+"isBelow"+str(block2)+str(shape2))
                    block2IsBelowBlock1 = model.NewBoolVar("blockShape"+str(block2)+str(shape2)+"isBelow"+str(block1)+str(shape1))

                    model.Add(xVarList[block1] + offset1[0] + offset1[2] <= xVarList[block2] + offset2[0]).OnlyEnforceIf(block1IsLeftOfBlock2)
                    model.Add(xVarList[block2] + offset2[0] + offset2[2] <= xVarList[block1] + offset1[0]).OnlyEnforceIf(block2IsLeftOfBlock1)
                    model.Add(yVarList[block1] + offset1[1] + offset1[3] <= yVarList[block2] + offset2[1]).OnlyEnforceIf(block1IsBelowBlock2)
                    model.Add(yVarList[block2] + offset2[1] + offset2[3] <= yVarList[block1] + offset1[1]).OnlyEnforceIf(block2IsBelowBlock1)

                    model.AddBoolOr([block1IsLeftOfBlock2, block2IsLeftOfBlock1, block1IsBelowBlock2, block2IsBelowBlock1])


    model.Minimize(lengthUsedVar)

    solver = cp_model.CpSolver()
    status = solver.Solve(model)

    if status == cp_model.INFEASIBLE:
        print("Infeasible model")

    if status == cp_model.OPTIMAL:
        print('Minimum length: %i' % solver.ObjectiveValue(), "\n")
        print("x =", [solver.Value(var) for var in xVarList])
        print("y =", [solver.Value(var) for var in yVarList])
    print("----------------------------")



#----------------------------------------------------------------------------------------------------------------------#




def rectilinearPackingWithRotationsProblem():
    model = cp_model.CpModel()

    # xoffset,yoffset,xsize,ysize
    offsetList = [(1,0,2,5), (0,1,5,2), (0,3,1,1), (1,0,2,5), (3,4,1,1), (3,1,1,1), (0,0,1,1), (3,3,1,1),
                  (4,0,1,1), (0,0,4,4), (1,4,3,1), (4,2,2,2), (0,2,4,4), (4,2,1,3), (2,0,2,2), (1,0,4,4),
                  (0,1,1,3), (1,4,2,2), (2,1,4,4), (2,0,3,1), (0,1,2,2), (1,5,1,2), (1,0,3,5), (4,1,1,4),
                  (5,3,2,1), (0,1,5,3), (1,0,4,1), (3,0,3,1), (0,1,2,1), (2,1,5,3), (2,4,4,1), (4,3,1,3),
                  (3,0,1,2), (1,2,3,5), (0,2,1,4), (0,0,1,3), (1,0,3,4), (0,3,3,1), (0,0,4,3), (1,0,3,1),
                  (0,1,4,3), (3,1,1,3), (0,0,3,4), (0,0,5,4), (0,0,4,5)]
    shapeList = [[(0,1,2), (2,3,4), (5,3,6), (7,1,8)],
                 [(9,10,11), (12,13,14), (15,16,17), (18,19,20)],
                 [(16,21,22,23), (10,24,25,26), (27,28,29,30), (31,32,33,34)],
                 [(43,), (44,)],
                 [(35,36), (37,38), (39,40), (41,42)]]
    numOfBlocks = len(shapeList)
    widthOfRiver = 9
    maxLength = 16


    xVarList = []
    yVarList = []
    for block in range(numOfBlocks):
        xVarList.append(model.NewIntVar(0, maxLength, "xPositionForBlock"+str(block)))
        yVarList.append(model.NewIntVar(0, widthOfRiver, "yPositionForBlock" + str(block)))

    lengthUsedVar = model.NewIntVar(0, maxLength, "lengthUsed")


    selectRotationVarList, selectRotationIndexDict, counter = [], {}, 0
    for block in range(numOfBlocks):
        varList = []
        for rotation in range(len(shapeList[block])):
            selectRotationVarList.append(model.NewBoolVar("block"+str(block)+"rotation"+str(rotation)))
            selectRotationIndexDict[(block, rotation)] = counter
            counter += 1
            varList.append(selectRotationVarList[-1])
        model.Add(sum(varList) == 1) # one rotation must be selected


    for block in range(numOfBlocks):
        for rotation, shapes in enumerate(shapeList[block]):
            for shape in shapes:
                offset = offsetList[shape]
                selectRotationIdx = selectRotationIndexDict[(block, rotation)]
                model.Add(xVarList[block] + offset[0] + offset[2] <= lengthUsedVar).OnlyEnforceIf(selectRotationVarList[selectRotationIdx])
                model.Add(yVarList[block] + offset[1] + offset[3] <= widthOfRiver).OnlyEnforceIf(selectRotationVarList[selectRotationIdx])


    for block1 in range(numOfBlocks-1):
        for rotation1, shapes1 in enumerate(shapeList[block1]):
            for shape1 in shapes1:
                offset1 = offsetList[shape1]

                for block2 in range(block1+1, numOfBlocks):
                    for rotation2,shapes2 in enumerate(shapeList[block2]):
                        for shape2 in shapes2:
                            offset2 = offsetList[shape2]

                            block1IsLeftOfBlock2 = model.NewBoolVar("blockShape"+str(block1)+str(shape1)+"isLeftTo"+str(block2)+str(shape2))
                            block2IsLeftOfBlock1 = model.NewBoolVar("blockShape"+str(block2)+str(shape2)+"isLeftTo"+str(block1)+str(shape1))
                            block1IsBelowBlock2 = model.NewBoolVar("blockShape"+str(block1)+str(shape1)+"isBelow"+str(block2)+str(shape2))
                            block2IsBelowBlock1 = model.NewBoolVar("blockShape"+str(block2)+str(shape2)+"isBelow"+str(block1)+str(shape1))
                            bothRotationsSelected = model.NewBoolVar("block1"+str(block1)+"rotation1"+str(rotation1)+"block2"+str(block2)+"rotation2"+str(rotation2)+"isSelected")

                            model.Add(xVarList[block1] + offset1[0] + offset1[2] <= xVarList[block2] + offset2[0]).OnlyEnforceIf(block1IsLeftOfBlock2)
                            model.Add(xVarList[block2] + offset2[0] + offset2[2] <= xVarList[block1] + offset1[0]).OnlyEnforceIf(block2IsLeftOfBlock1)
                            model.Add(yVarList[block1] + offset1[1] + offset1[3] <= yVarList[block2] + offset2[1]).OnlyEnforceIf(block1IsBelowBlock2)
                            model.Add(yVarList[block2] + offset2[1] + offset2[3] <= yVarList[block1] + offset1[1]).OnlyEnforceIf(block2IsBelowBlock1)

                            selectRotation1Idx, selectRotation2Idx = selectRotationIndexDict[(block1, rotation1)], selectRotationIndexDict[(block2, rotation2)]
                            model.AddMultiplicationEquality(bothRotationsSelected, [selectRotationVarList[selectRotation1Idx], selectRotationVarList[selectRotation2Idx]])
                            model.Add(block1IsLeftOfBlock2 + block2IsLeftOfBlock1 + block1IsBelowBlock2 + block2IsBelowBlock1 >= bothRotationsSelected)


    model.Minimize(lengthUsedVar)

    solver = cp_model.CpSolver()
    status = solver.Solve(model)

    if status == cp_model.INFEASIBLE:
        print("Infeasible model")

    if status == cp_model.OPTIMAL:
        print('Minimum length: %i' % solver.ObjectiveValue(), "\n")
        print(solver.Value(lengthUsedVar))
        print("x =", [solver.Value(var) for var in xVarList])
        print("y =", [solver.Value(var) for var in yVarList])
    print("----------------------------")



#----------------------------------------------------------------------------------------------------------------------#





if __name__ == "__main__":

    squarePackingProblem()
    rectilinearPackingProblem()
    rectilinearPackingWithRotationsProblem()