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)