46、动态规划-机器更新问题

import random as rd
from pyscipopt import Model, quicksum

def C_generate(n, dC):

    C = [0 for i in range(n)]

    for i in range(n-1):
        C[i+1] = C[i] + int(rd.random() * dC)

    return C

def O_generate(n, dO):

    O = [0 for i in range(n)]

    for i in range(n-1):
        O[i+1] = O[i] + int(rd.random() * dO)

    return O

def I_generate(n, I0, dT):

    I = [I0 for i in range(n)]

    for i in range(n - 1):
        I[i + 1] = I[i] - int(rd.random() * dT)

    return I

def planing(C,O,I, time_limit):

    M = 10000

    n = len(C)

    dO = [O[i] - O[i - 1] if i > 0 else O[i] for i in range(n)]
    dI = [I[i] - I[i - 1] if i > 0 else I[i] for i in range(n)]
    print("\ndO:\n", dO)
    print("\ndI:\n", dI)

    model = Model("Machine_update")

    #A[i][j]代表第i年使用的机器役龄超过j年,U[i]==1代表第i年更新机器
    A = [[model.addVar(vtype="B", name="A[%s,%s]" % (i, j)) for j in range(n)] for i in range(n)]
    U = [model.addVar(vtype="B", name="U[%s]" % (i)) for i in range(n)]

    # 以总成本最小为目标
    model.setObjective(quicksum((dI[j]-dO[j])*A[i][j] for i in range(n) for j in range(n)) - quicksum(C[i]*U[i] for i in range(n)), "maximize")

    #服役j+1年得先服役j年
    for i in range(n):
        for j in range(n-1):
            model.addCons(A[i][j] - A[i][j+1] >= 0)

    #更新时的机器服役0年
    for i in range(n):
        model.addCons(A[i][0] + U[i] == 1)

    #最开始机器是新的
    model.addCons(A[0][0] == 0)

    #机器要么换新的，要不年龄长一岁
    for i in range(n-1):
        model.addCons(quicksum(A[i + 1][j] for j in range(n)) - quicksum(A[i][j] for j in range(n)) - 1 + M * U[i + 1] >= 0)
        model.addCons(quicksum(A[i + 1][j] for j in range(n)) - quicksum(A[i][j] for j in range(n)) - 1 - M * U[i + 1] <= 0)
        model.addCons(quicksum(A[i + 1][j] for j in range(n)) - M * (1-U[i + 1]) <= 0)
        model.addCons(quicksum(A[i + 1][j] for j in range(n)) + M * (1-U[i + 1]) >= 0)

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

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

    return U1,A1

if __name__ == "__main__":

    # 年数
    n = 5

    # 每个年初更新机器产生的费用（假设最开始时机器是新的）
    C = C_generate(n, 4)
    print("\nC:\n", C)

    # 运行费用
    O = O_generate(n, 1547)
    print("\nO:\n", O)

    #机器产生的价值
    I = I_generate(n, 20, 5)
    print("\nI:\n", I)

    #更新策略生成
    U,A = planing(C,O,I,100)
    print("\nU:\n", U)
    print("\nA:\n", A)