pulp_with_initial_solution.py

from pulp import *


def solve(initial_solution_file_path: str=None):
    left = list(range(10))
    right = list(range(10))
    route = [(i, j) for i in left for j in right]

    supply = [8, 6, 7, 10, 11, 13, 6, 13, 14, 12]
    demand = [11, 11, 10, 7, 10, 11, 13, 9, 9, 9]
    assert(sum(supply) == sum(demand))

    # 枝(i, j)を使う定数コスト
    costs = [[1] * len(right) for _ in range(len(left))]
    costs = makeDict([left, right], costs, 0)

    M = [[0] * len(right) for i in range(len(left))]
    for i in left:
        for j in right:
            M[i][j] = min(supply[i], demand[j])
    M = makeDict([left, right], M, 0)

    x = LpVariable.dicts(name="flow", indexs=route, lowBound=0, cat=LpInteger)  # 枝(i,j)に流れる量
    y = LpVariable.dicts(name="use", indexs=route, cat=LpBinary)    # 枝(i,j)を使うかどうか

    prob = LpProblem("PC-FCTP", LpMinimize)
    prob += lpSum([costs[r[0]][r[1]] * y[r] for r in route])

    for i in left:
        prob += lpSum([x[(i, j)] for j in right]) == supply[i]

    for j in right:
        prob += lpSum([x[(i, j)] for i in left]) == demand[j]

    for r in route:
        prob += x[r] <= M[r[0]][r[1]] * y[r]

    options = None
    if initial_solution_file_path is not None:
        options = ["mips", initial_solution_file_path]

    prob.writeLP("problem.lp")
    prob.solve(PULP_CBC_CMD(options=options, msg=1))

    print("Status:{0}".format(LpStatus[prob.status]))
    print("Total Cost:{0}".format(value(prob.objective)))
    print("Assignment:")
    for r in route:
        if y[r].value() == 1:
            print("Left:{0} -> Right:{1}({2})".format(r[0], r[1], x[r].value()))


def make_initial_solution(initial_solution_file_path):
    # 初期解
    d = {(0, 4): 8,
         (1, 1): 6,
         (2, 3): 7,
         (3, 2): 10,
         (4, 0): 11,
         (5, 6): 13,
         (6, 8): 6,
         (7, 4): 2,
         (7, 5): 11,
         (8, 1): 5,
         (8, 7): 9,
         (9, 8): 3,
         (9, 9): 9}

    left = list(range(10))
    right = list(range(10))
    route = [(i, j) for i in left for j in right]

    supply = [8, 6, 7, 10, 11, 13, 6, 13, 14, 12]
    demand = [11, 11, 10, 7, 10, 11, 13, 9, 9, 9]
    assert(sum(supply) == sum(demand))

    # 枝(i, j)を使う定数コスト
    costs = [[1] * len(right) for _ in range(len(left))]
    costs = makeDict([left, right], costs, 0)

    M = [[0] * len(right) for i in range(len(left))]
    for i in left:
        for j in right:
            M[i][j] = min(supply[i], demand[j])
    M = makeDict([left, right], M, 0)

    x = LpVariable.dicts(name="flow", indexs=route, lowBound=0, cat=LpInteger)  # 枝(i,j)に流れる量
    y = LpVariable.dicts(name="use", indexs=route, cat=LpBinary)                # 枝(i,j)を使うかどうか

    prob = LpProblem("PC-FCTP", LpMinimize)
    prob += lpSum([costs[r[0]][r[1]] * y[r] for r in route])

    for i in left:
        prob += lpSum([x[(i, j)] for j in right]) == supply[i]

    for j in right:
        prob += lpSum([x[(i, j)] for i in left]) == demand[j]

    for r in route:
        prob += x[r] <= M[r[0]][r[1]] * y[r]

    # 初期解
    for r in route:
        if r in d:
            v = d[r]
            prob += x[r] == v
            prob += y[r] == 1
        else:
            prob += x[r] == 0
            prob += y[r] == 0

    prob.solve(PULP_CBC_CMD(options=["solve", "solution", initial_solution_file_path]))
    print("Status:{0}".format(LpStatus[prob.status]))
    print("Total Cost:{0}".format(value(prob.objective)))


def main():
    import time
    solution_file = "initialSolution.sol"

    print("make initial solution")
    make_initial_solution(initial_solution_file_path=solution_file)
    start = time.time()

    print("solve")
    solve(initial_solution_file_path=solution_file)
    print("{0}sec".format(int(time.time() - start)))


if __name__ == '__main__':
    main()