42、动态规划-资源分配问题

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

#资源分配问题，例如n台机器分给m个工厂

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 G_generate(m, n):

    G = [[0 for j in range(n+1)] for i in range(m)]

    for i in range(m):
        for j in range(1,n+1):
            G[i][j] = G[i][j-1] + int(rd.random() * 5)

    return G

def IP(dG, time_limit):

    m = len(dG)
    n = len(dG[0])-1

    model = Model("Resource_allocation_IP")
    A = [[model.addVar(vtype="B", name="A[%s,%s]" % (i, j)) for j in range(n+1)] for i in range(m)]

    # 以总成本最小为目标
    model.setObjective(quicksum(dG[i][j]*A[i][j] for i in range(m) for j in range(n+1)), "maximize")

    # A[i][j] >= A[i][j+..]
    for i in range(m):
        for j in range(n):
            for k in range(j+1, n+1):
                model.addCons(A[i][j] - A[i][k] >= 0)

    model.addCons(quicksum(A[i][j] for i in range(m) for j in range(1,n+1)) == n)

    # 设置求解时间
    model.setRealParam("limits/time", time_limit)
    model.optimize()
    print("\ngap:",model.getGap())

    # 拿结果
    A1 = [[round(model.getVal(A[i][j])) for j in range(n + 1)] for i in range(m)]

    return A1

def DP(dG):

    m = len(dG)
    n = len(dG[0])-1

    # 局部结果，A[i][j]表示j个机器分给第i个工厂的情况下，
    A = [[0 if i > 0 else sum(dG[i][:j+1]) for j in range(n+1)] for i in range(m)]
    A_detail = [[[0 if i > 0 or k > 0 else j for k in range(m)] for j in range(n+1)] for i in range(m)]

    for i in range(1,m):
        gi = [sum(dG[i][:j+1]) for j in range(n+1)]
        for j in range(1, n+1):
            idex = 0
            for k in range(1, j+1):
                if A[i-1][idex] + gi[j - idex] < A[i-1][k] + gi[j - k]:
                    idex = k
            A[i][j] = A[i-1][idex] + gi[j - idex]
            A_detail[i][j] =  [A_detail[i-1][idex][k] if k != i else j - idex for k in range(m)]

    R = A_detail[m-1][n]

    A1 = [[1 if j <= R[i] else 0 for j in range(n+1)] for i in range(m)]

    return A1

def allocation(dG, time_limit = 100):

    A_IP = IP(dG, time_limit)
    A_DP = DP(dG)

    return A_IP, A_DP

if __name__ == "__main__":

    #工厂数、机器数
    m,n = 3,5

    #每个工厂有不同机器数目能产生的盈利
    G = G_generate(m, n)
    LL_print("G", G)

    #每台机器的边际效用
    dG = [[0 if j == 0 else G[i][j] - G[i][j - 1] for j in range(n + 1)] for i in range(m)]
    LL_print("dG", dG)

    # 求解，得到每个工厂分配的机器数
    A_IP, A_DP = allocation(dG, 100)
    LL_print("A_IP", A_IP)
    LL_print("A_DP", A_DP)

    # 每个工厂分配的机器数汇总
    R_IP = [sum(A_IP[i][1:]) for i in range(m)]
    R_DP = [sum(A_DP[i][1:]) for i in range(m)]
    print("\nR_IP:\n", R_IP)
    print("\nR_DP:\n", R_DP)

    # 价值
    g_IP = sum([sum(dG[i][:(R_IP[i] + 1)]) for i in range(m)])
    g_DP = sum([sum(dG[i][:(R_DP[i] + 1)]) for i in range(m)])
    print("\ng_IP:\n", g_IP)
    print("\ng_DP:\n", g_DP)