动态规划求解带移动平均约束的多周期产品最大销售额调度问题

问题说明 #

我们需要在T个周期内销售总量为N的产品，拆分各周期计划销量为n₀,n₁,⋯,n_{T−1}，满足∑nᵢ=N，目标是最大化总实际销售额∑Sᵢ。
实际销售额Sᵢ的计算依赖当期移动平均M和当期计划销量nᵢ，给定参数α=0.001、π=0.5，计算规则如下：

• 初始化M=0，按i=0,1,…,T−1顺序迭代计算
• 当期移动平均更新规则：Mᵢ = ⌈0.5*(上期M值 + nᵢ)⌉ （⌈⌉代表向上取整）
• 当期实际销售额计算规则：Sᵢ = ⌈(1−α*Mᵢ^π)*nᵢ⌉

正向计算参考代码 #

import math
import numpy as np

M = 0
T = 4
N = 10000
alpha = 1e-3
pi = 0.5
S = np.zeros(T,dtype='i')
n  = np.array([5000,1000,2000,2000])
print(n)
total = 0
for i in range(T):
    M = math.ceil(0.5*(M + n[i]))
    S[i] = math.ceil((1 - alpha*M**pi)*n[i])
    total += S[i]
    print('at time %d, M = %d and we sell %d products' %(i,M,S[i]))
print('total sold =', total)

原实现存在的问题 #

原思路采用三维状态存储倒推结果的逻辑是成立的，问题出在以下几个细节：

1. 销售额计算逻辑顺序错误：按照规则需要先更新当期移动平均M，再用新的M计算当期销售额，原代码直接用期初的m值计算当期销售额，和规则不符
2. 状态越界风险：原代码假设移动平均m的上界为当前剩余库存n，当更新后的M值大于下一周期的剩余库存n-plan时，访问下一周期状态会出现数组越界，正确的上界应该是总库存N（M最大值不可能超过总库存N）
3. 边界状态计算错误：最后一个周期的销售额计算同样需要先更新M，再算销售额，原代码直接用输入的m值计算

修复后的代码 #

import math
import numpy as np

def DP_opt(N, T, alpha, pi, dp):
    # 边界：最后一个周期，把剩余所有库存都卖掉
    for n_remain in range(0, N+1):
        for m_start in range(0, N+1):
            # 先更新M，再算销售额
            m_new = math.ceil(0.5 * (m_start + n_remain))
            dp[T-1][n_remain][m_start] = math.ceil((1 - alpha * (m_new ** pi)) * n_remain)
    
    # 倒推前面的周期
    for t in range(T-2, -1, -1):
        print(f"处理周期{t}")
        for n_remain in range(0, N+1):
            for m_start in range(0, N+1):
                best = 0
                # 当期计划销量plan的取值范围是0到剩余库存n_remain
                for plan in range(0, n_remain + 1):
                    # 先算更新后的M
                    m_new = math.ceil(0.5 * (m_start + plan))
                    # 用新的M算当期销售额
                    sale_now = math.ceil((1 - alpha * (m_new ** pi)) * plan)
                    # 下一周期的剩余库存是n_remain - plan，期初M是m_new
                    sale_later = dp[t+1][n_remain - plan][m_new]
                    total_candidate = sale_now + sale_later
                    if total_candidate > best:
                        best = total_candidate
                dp[t][n_remain][m_start] = best
    
    # 初始状态：第0个周期，剩余库存N，期初M=0
    print("最大总销售额：", dp[0][N][0])
    return dp

# 调用示例
if __name__ == "__main__":
    N = 100
    T = 5
    pi = .5
    alpha = 1e-3
    # dp维度：周期数 * 剩余库存+1 * 移动平均最大值+1
    dp = np.zeros((T, N+1, N+1), dtype=np.int32)
    DP_opt(N, T, alpha, pi, dp)

内容的提问来源于stack exchange，提问作者Cindy Philip